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SECTION  I 


INTRODUCTION 

With  recent  developments  in  target  hardening  technology*  damage  to  the 
casings  of  impacting  bombs  prior  to  initiation  of  the  explosive  charge  by  the 
fuse  has  been  occurring  with  increasing  frequency.  The  work  done  here  was 
motivated  by  the  desire  to  be  able  to  predict  whether  or  not  detonation  will 
occur  in  warheads  whose  casings  have  been  fractured.  The  existence  of  a  hole 
in  the  casing  which  confines  the  explosive  can  be  expected  to  have  a 
considerable  effect  on  the  formation  of  a  detonation  wave  in  the  explosive. 

The  purpose  of  this  work  is  to  present  a  model  which  will  determine  the  extent 
of  this  effect. 

An  on-going  effort  has  been  present  at  the  University  of  Illinois  at 
Urbana-Charapaign  for  the  past  several  years,  under  the  direction  of  Professor 
Herman  Krier,  to  develop  a  model  which  accurately  describes  the  fluid 
mechanics  that  result  from  flame  spreading  through  a  fragmented  (and  thus 
porous)  high-energy  solid  propellant  or  explosive  media.  The  thrust  of  this 
effort  has  been  to:  (a)  formulate  the  equations  of  motion  governing  the  fluid 
dynamics  of  a  combusting  propellant  bed  and,  (b)  develop  a  stable  numerical 
scheme  using  these  equations  that  will  predict  the  transition  from 
deflagration  to  detonation  (DDT)  in  a  confined  propellant  bed.  Previous 
studies  (including  those  of  Van  Tassel  and  Krier  ill,  Krier  and  Gokhale  [2], 
Krier,  Rajan,  and  Van  Tassel  (3),  Dimitsein  [4],  Krier,  Dimitstein,  and 
Gokhale  [5],  Krier,  Gokhale,  and  Hughes  |6l,  Krier  and  Kezerle  [7],  and 
Butler,  Krier,  and  Lembeck  181)  have  laid  a  strong  foundation  towards  the 
achievement  of  the  goal.  Models  were  derived  to  treat  the  unsteady,  two- 


1 


phase,  separated  flow  conservation  equations  in  one  space  dimension  and  these 
were  incorporated  into  a  computer  code  which  yielded  stable  solutions  for  a 
wide  range  of  input  parameters.  In  its  present  form,  the  existing  code  is 
capable  of  showing  the  rapid  building  up  of  a  pressure  shock  front  from  an 
initially  quiescent,  but  locally  ignited  bed  of  packed,  granulated  particles 
and  predicting  the  formation  of  a  steady-state  detonation  wave,  provided  that 
the  motion  of  the  gas  and  particle  phases  is  in  one  space  dimension  only. 

(This  has  been  shown  experimentally  to  be  the  case  for  flame  spreading  through 
a  granulated  propellant  or  explosive  material  which  is  completely  enclosed  in 
a  solid  container  (Reference  5).)  Table  I  is  a  list  which  compares  the 
detonation  conditions  predicted  by  this  one-dimensional  ODT  code  for  various 
values  of  the  initial  particle  density  and  porosity  with  the  results  predicted 
by  a  steady-state  thermo-chemical  code  (Reference  10);  it  illustrates  the 
high  degree  of  accuracy  using  these  methods. 


TABLE  1.  COMPARISON  OF  DATA  FROM  1-D  CODE  AND  TIGER  CODE 
(oQ  AND  *0  ARE  INPUT  PARAMETERS) 


p0(g/cc) 

♦o 

PCJ(GPa) 

TcjCM 

pCJ(g/cc) 

DDT  -  Code 

1.20 

0.368 

14.38 

4201 

1.64 

TIGER 

1.20 

0.368 

14.92 

4337 

1.65 

ODT  -  Code 

1.30 

0.316 

16.90 

4289 

1.76 

TIGER 

1.30 

0.316 

17.26 

4304 

1.78 

DOT  -  Code 

1.33 

0.300 

18.11 

4406 

1.80 

TIGER 

1.33 

0.300 

18.00 

4300 

1.81 

DDT  -  Code 

1.40 

0.263 

19.64 

4393 

1.87 

TIGER 

1.40 

0.263 

19.60 

4280 

1.89 
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In  its  current  state,  the  model  can  predict  DDT  in  fractured  high 
explosives  contained  within  undamaged  bomb  casings,  but  does  not  contain  the 
necessary  considerations  for  modeling  the  accelerated  flame  spreading  in 
realistic  2-D  and  3-D  bomb  casing  configurations,  since  one  of  the  assumptions 
in  its  formulation  was  that  the  flow  was  strictly  one  dimensional.  The 
addition  of  cracks  or  holes  in  the  container  walls  (partial  confinement)  would 
cause  the  flow  to  become  multi-dimensional.  The  work  presented  in  this  report 
undertakes  to  develop  a  model  of  a  combusting  propellant  bed  which  is  not 
always  totally  confined,  that  is,  in  which  some  of  the  gas  generated  by  the 
burning  fragments  is  allowed  to  escape  (along  with,  perhaps,  some  of  the 
entrained  particles). 

Two  approaches  to  solving  this  problem  are  discussed:  (a)  modifying  the 
conse  /ation  equations  in  the  one-dimensional  model  to  include  loss  terms 
which  simulate  the  mass,  momentum,  and  energy  escaping  through  the  holes  in 
the  casing,  and  (b)  developing  a  new,  unsteady,  three-dimensional  flow  model 
with  a  new  numerical  scheme,  which  incorporates  the  hole  in  the  casing  as  a 
boundary  condition  to  the  flow. 

Approach  (a),  which  is  discussed  in  detail  in  Section  II,  utilizes 
pseudo-sink  terms  in  the  conservation  equations  which  simulate  the  effects  of 
mass  loss  on  the  formation  of  a  pressure  wave  in  the  burning  particle  bed. 

The  approach  was  taken  as  a  first  approximation  of  the  effect  of  partial 
confinement  since  the  working  one-dimensional  model  was  already  available. 
Clearly,  since  the  mass  would  be  ejected  out  of  the  domain  in  a  direction 
transverse  to  the  direction  of  the  flame  spreading,  the  problem  becomes  multi¬ 
dimensional.  However,  even  though  the  quasi  two-dimensional  formulation  will 
not  yield  exact  results,  it  can  provide  meaningful  physical  insight  into  the 


behavior  of  partially  confined  accelerating  reactive  flows.  This  phase  of  the 
work  has  been  successfully  completed  and  published  in  a  proceedings  (Reference 
11).  The  results  are  both  reasonable  and  consistent. 

Approach  (b)  is  much  more  ambitious  in  its  scope.  It  requires  the  two- 
phase,  unsteady,  separated  flow  conservation  equations  to  be  expressed  in 
three  space  dimensions  and  a  new  stable  numerical  scheme  to  be  developed  which 
can  solve  these  very  non-linear,  coupled  time-dependent  equations.  Once  this 
has  been  done,  the  condition  of  partial  confinement  can  then  be  added  in  by 
simply  using  the  appropriate  boundary  condition.  However,  the  difficulty  of 
developing  a  stable  numerical  scheme  to  integrate  these  coupled,  hyperbolic 
partial  differential  conservation  equations  alone  provides  a  formidable 
task.  When  coupled  to  the  large  domain  of  the  problem  the  task  becomes  even 
more  difficult.  It  is  therefore  apparent  that  this  phase  of  the  work  cannot 
be  expected  to  have  been  seen  through  to  its  completion  in  the  short  year  that 
it  has  been  underway.  However,  important  strides  have  been  made  toward  that 
completion.  Section  III  of  this  report  presents  the  formulation  of  the 
unsteady,  three-dimensional,  two-phase  reactive  flow  model,  while  the 
numerical  scheme  is  discussed  in  Section  IV.  Results  currently  available  are 
presented  in  Section  V. 
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SECTION  II 


THE  PSEUDO  TWO-DIMENSIONAL  FORMULATION 


1.  INTRODUCTION 

Impact  forces  imparted  from  hardened  concrete  targets  to  aerial  bombs 
can,  in  some  cases,  cause  outer  case  failure  and  breakup  of  the  high  explosive 
filler.  The  detonation  wave  which  may  be  expected  to  form  if  the  bomb  casing 
were  undamaged  (total  confinement  of  the  explosive  media)  may  be  reduced  or 
even  quenched  completely  in  the  presence  of  the  damaged  casing  (partial 
confinement).  In  this  chapter  the  effect  of  partial  confinement  is 
investigated  by  adding  quasi  two-dimensional  sink  terms  to  the  equations  of 
mass,  momentum,  and  energy  that  describe  the  accelerating  flame  spreading 
through  a  granular  high  explosive.  These  sink  terms  are  formulated  by 
assuming  that  the  gas  escapes  from  the  hole  at  the  choked  velocity  for  the 
local  pressure  conditions  at  the  hole  location.  It  is  clear  from  the  results 
presented  that  such  partial  confinement  can  indeed  have  a  significant  effect 
on  the  formation  of  a  detonation  wave. 

2.  Assumptions 

The  important  assumptions  in  formulating  this  psuedo  one-dimensional 
model  are  listed  below. 

a.  The  analysis  considered  here  explicitly  assumes  that  at  time  t  =  0, 
the  warhead  has  already  impacted  with  the  surface,  and  the  explosive  filler 
has  already  fragmented,  causing  damage  to  the  casing  (see  Figure  1).  It  is 
also  assumed  that  at  t  =  QJ,  ignition  of  a  small  region  of  the  explosive  by  the 
fuse  had  already  occurred. 
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b.  The  deflagration  and  possible  transition  to  detonation  occurs  only 
in  one  direction  (x-space). 

c.  The  surface-to-volume  ratio  of  the  fragmented  particles  is 
sufficiently  high  that  they  may  be  treated  as  uniform  pseudo-spheres, 
typically  of  the  millimeter  or  sub-millimeter  diameter. 

d.  The  unsteady  two-phase,  separated  flow  analysis  previously  developed 
in  Reference  8  represents  the  basis  upon  which  this  model  was  built. 

e.  The  possibility  of  mass,  momentum,  and  energy  loss  is  treated  by 
psuedo-side  venting  from  cracks  of  prescribed  width,  length,  and  location. 

f.  The  decision  on  whether  the  fragmented  bed  of  explosive  will 
detonate  is  based  on  the  transient  reacting  flow  events  occurring  in  the  first 
10  to  20  cm  of  length.  A  run-up  length  of  greater  than  20  cm  requires  times 
well  beyond  established  experimental  DOT  events. 

g.  The  explosive  particles,  when  ignited,  burn  at  a  rate  f  =  aPn  where 
a  and  n  are  assumed  to  be  known  constants  for  typical  high  explosives  used  in 
warheads.  Ignition  is  assumed  to  have  occurred  when  a  prescribed  critical 
temperature  is  reached  by  the  particles  (see  Reference  8). 

h.  The  porosity  (gas  volume/total  volume)  of  the  fragmented  bed  is 
assumed  initially  to  be  uniform,  and  typically  in  the  range  of  0.2  to  0.3. 
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i.  The  product  gases  obey  a  non-ideal  equation  of  state  (see  Reference 

7). 

j.  Heat  transfer  from  the  gas  to  the  solid  is  by  rapid  convection  only; 
conduction  and  radiation  are  ignored. 

k.  The  gas  is  assumed  to  be  inviscid.  The  drag  between  the  gas  and 
particle  phase  is  included  as  a  source/sink  term  in  the  momentum  and  energy 
equations.  (This  is  equivalent  to  assuming  that  the  gas-gas  momentum  viscous 
losses  are  orders  of  magnitude  less  than  the  gas-particTe  viscous 
interaction.) 

l.  The  equations  are  expressed  as  averaged- laminar  flow  properties  in 
that  the  turbulence  due  to  the  two-phase  nature  of  the  problem  has  been 
averaged  out. 

3.  THE  FLUID  MECHANICS  MODEL 

In  the  analysis  of  two-phase  reactive  flow,  one  must  describe  the 
conservation  of  mass,  momentum,  and  energy  throughout  the  domain  for  both  the 
solid  particle  phase  and  the  gaseous  products  phase.  In  separated  flow 
analysis,  each  phase  is  assumed  to  be  a  continuum,  and  the  important 
quantities  of  mass,  momentum,  and  energy  are  therefore  conserved  separately  In 

each  phase.  The  governing  differential  equations  are  presented  below, 

/ 

including  the  modifications  made  in  this  work  to  account  for  the  ldsses 
occurring  through  the  hole  in  the  casing.  For  a  detailed  derivation  of  the 
basic  equations,  the  reader  is  referred  to  Reference  7. 
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In  the  energy  equation, 


and 


:pt 


-  C  T 
Lvp'p 


+  0.5  U, 


The  subscripts  g  and  p  denote  the  gas  and  particle  phase  respectively.  The 
quantities  p^  and  p2  equations  (1)  -  (6)  refer  to  the  bulk  densities  for 
each  phase  and  are  defined  as 


Pl  =  Og$ 

°2  ~ 

The  porosity  $  is  defined  as  the  ratio  of  instantaneous  gas  volume  to  the 
mixture  volume.  Hence,  the  solids  fraction  is  l-<>.  The  8g  and  ep  terms  in 
Equations  (1)  -  (6)  are  the  sink  terms  which  account  for  the  losses  through 
the  cracked  casing.  They  are  derived  in  detail  in  the  next  section. 

In  addition  to  the  six  conservation  equations,  several  constitutive 
relations  are  necessary  for  closure.  The  relations  used  in  this  work  are 
presented  in  detail  in  Appendix  A,  along  with  some  discussion  on  their 
validity  under  the  extreme  conditions  imposed  by  this  problem.  The  necessary 
equations  are: 

(A-l)  An  equation  of  state  for  the  gaseous  products  phase,  Pg  =  Pg(og,  Eg). 

(A-Z)  An  equation  of  state  for  the  solid  particle  phase,  pp  =  op(Pp). 

(A-3)  A  relation  for  predicting  the  instantaneous  porosity,  <j>. 

(A-4)  A  relation  for  the  interphase  heat  transfer,  Q  . 

(A-5)  A  gas  particle  interphase  drag  relation,  0. 

(A-6)  An  equation  determining  the  gas  generation  by  the  burning  particles,  r 
(A-7)  An  ignition  criteria. 
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4.  DERIVATION  OF  THE  MASS  LOSS  TERM 

The  condition  of  partial  confinement  requires  that  extra  terms  be  added 
to  the  conservation  equations  in  Reference  7.  Figure  2  shows  the  fragmented 
bed  being  modeled.  The  fourth  volume  cell  from  the  left  shows  the  crack  in 
the  casing,  exposing  the  granulated  propellant  inside.  Superimposed  on  the 
front  face  of  the  illustration  are  the  pressure  versus  distance  profiles 
showing  the  expected  pressure  drop  at  the  hole. 

Consider  now  the  single  control  volume  with  a  hole  in  it  shown  in  Figure 
3.  Pressure  gradients  aue  to  combustion  are  causing  gas  flow  in  the  x- 
direction,  as  well  as  causing  some  of  the  mass  to  flow  out  of  the  hole.  The 
rate  of  mass  flow  of  the  gas  out  of  the  hole,  per  unit  volume,  is  given  by: 


s 


g 


\  =  °gUAh 


(7) 

(8) 


where  og  =  gas  density 

U  =  velocity  of  gas  leaving  the  control  volume 
Ah  =  area  of  the  hole 

Due  to  the  high  gas  pressures  present  in  the  bed,  we  can  assume  that  the 
flow  of  the  gas  out  of  the  hole  into  the  atmosphere  will  always  be  choked. 
Hence,  U  will  always  be  the  local  speed  of  sound: 


(9) 
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A  detailed  derivation  of  the  speed  of  sound  equation  for  the  non-ideal  gas 
used  in  this  model  can  he  found  in  Appendix  B. 


The  area  of  the  hole  is  defined  as: 


=  ax 'ax* 


(10) 


Note  that  *  is  included  in  the  definition  of  Ah  since,  in  the  separated  flow 
analysis,  only  the  area  of  the  hole  occupied  by  the  gas  is  considered. 

The  volume  of  the  control  cell  is  given  by  Vcy  =  ax^.  Substituting  back 
into  the  expression  for  gg,  one  obtains 


oqU$ax  flx'* 
_  =  _a  - 

cv  ax 


V 


ax 


(H) 

(12) 


where  CQ  is  a  prescribed  constant  which  defines  the  width  of  the  cack. 

If  one  were  to  consider  entrainment  of  some  of  the  particles  of  unburned 
solid  propellant  in  the  gas  which  is  escaping  out  of  the  hole,  the  mass  loss 
for  the  solid  phase  may  be  written  as: 


°pUpAh 
6n  =  -J- 

H  ax 


P,U  0 

3  =  -Ls_o  c 
p  ax  o 


(13) 


where  D0  is  a  prescribed  drag  factor  between  the  entrained  particles  and  the 
escaping  gas,  0  <  0Q  <  1. 


The  momentum  losses  through  the  hole  then  become  BgUg  and  epUp  for  the 
gas  and  particle  phase,  respectively.  Note  that  the  velocity  in  the  /- 
direction  is  used  in  these  terms  rather  than  the  sound  speed.  This  is  done 
because  of  the  one-dimensional  formulation  of  the  equations  of  motion.  The 
quantity  of  interest  is  the  amount  of  momentum  in  the  x-direction  which  the 
ejected  mass  removes  from  the  domain. 

Likewise,  the  sink  terms  in  the  energy  equations  are  BgEgt  and  spEpt  for 
the  gas  and  particle  phase,  respectively.  Note  that  the  work  term  found  in 
the  flux  terms  of  Equations  (5)  and  (6)  are  not  included  in  the  mass  loss  term 
because  there  is  no  work  done  in  the  x-direction  in  removing  the  mass  from  the 
domain. 


5.  NUMERICAL  SOLUTION  OF  THE  CONSERVATION  EQUATIONS 

The  highly  non-linear  partial  differential  equations  (1)  -  (6)  do  not 
admit  any  analytical  solutions.  In  order  to  obtain  solutions,  the  time- 
dependent  equations  were  discretized  over  the  domain  and  solved  numerically 

r 

via  the  method  of  lines.  In  this  approach,  the  flux  terms  on  the  right-hand 
side  of  Equations  (I)  -  (6)  are  evaluated  as  centered  finite  differences, 
along  with  any  source/sink  terms.  The.  remaining  differential  equations  are 
then  treated  as  ordinary  differential  equations  and  solved  via  an  O.D.E. 
solver.  For  example.  Equation  (1): 


3P, 

at 


»(p.V 

3X 


+  r 


(14) 


would  be  solved  by  first  evaluating  the  terms  on  the  right-hand  side  of  the 
equation  at  the  current  time  level  to  get,  say,  G(t).  This  result  is  then  put 
back  into  Equation  (1)  to  give: 
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(15) 


This  is  an  ordinary  differential  equation  with  as  the  dependent  variable 
and  t  as  the  only  independent  variable.  Now  an  O.D.E.  solver  such  as  Runge- 
Kutta  may  be  employed  to  predict  the  new  value  of  p^.  In  this  work  a  variable 
time  step  Runge-Kutta-Fehlberg  fourth-fifth  order  algorithm  was  used.  (For 
more  information  on  the  method  of  lines,  the  reader  is  referred  to  Reference 
12.) 

6.  COMPUTEO  RESULTS 

The  model  described  in  the  previous  sections  is  first  applied  to  a 
baseline  case,  one  in  which  there  is  no  hole  in  the  confining  wall  (total 
confinement,  CQ  =  0.0).  The  results  for  this  case  are  then  compared  with  the 
output  for  cases  where  holes  of  various  length,  width,  and  location  have  been 
introduced  into  the  wall.  The  important  input  parameters  for  these  cases  are 
listed  in  Table  2.  The  output  presented  here  is  primarily  in  the  form  of 
pressure  profiles  and  plots  of  the  locus  of  the  flame  front.  These  profiles 
are  presented  in  preference  to  others  because  pressure  and  flame  front 
velocity  are  the  two  most  critical  parameters  in  determining  whether  or  not  a 
detonation  will  occur. 

Figure  4  shows  the  pressure  versus  distance  profiles  for  various  times 
for  the  baseline  (no  hole)  case.  The  plot  shows  the  shock  building  and 
advancing  through  the  bed  as  time  progresses.  At  x  =  9.0  cm  the  peak  pressure 
is  almost  14.0  GPa,  which  is  just  a  little  less  than  the  Chapman-Jouget  (CJ) 
detonation  pressure  for  an  explosive  with  this  bulk  density  (see  Reference  8). 
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TABLE  2.  INPUT  PARAMETERS  FOR  1-0  MODEL 


po 


Bed  Length,  L 
Particle  Radius, 

Porosity, 

Pressure,  Pao 
Gas  Temperature,  T 
Particle  Temperature,  T 
Ignition  Temperature,  T1a 
Chemical  Energy  of  Solia,  E 


po 


chem 


Gas  Constant,  R' 

Prandtl  Number  of  Gas,  Pr 
Constant  Volume  Specific  Heat,  Cv 
Interphase  Viscosity,  u 
Burning  Rate  Index,  n 
Coefficient  in  Burning  Rate,  a 
Co-efficient  in  Equation  of  State, 
Node  Size,  ax 


10  cm 
200  ym 
0.3  _ 

1  x  105  GPa 
300  °K 
300  °K 
306  SK 
5.74  MJ/kg 
296.8  J/kg  °K 
0.70 

1500  J/kg  °K 
1.80  x  10“3  Ns/rrr 
1.0 
0.001 
4.0 
1  mm 


Figure  5  shows  the  effect  of  having  placed  a  hole  4  mm  in  length  with  an 
opening  ratio  (CQ)  of  0.5  at  a  distance  of  2.5  cm  from  the  initiated  end  of 
the  bed.  The  location  of  the  hole  on  the  plot  is  evident  from  the  severe 
indentation  in  the  pressure  profile.  It  is  interesting  to  compare  the  profile 
at  t  =  68  ysec  for  Figures  4  and  5.  As  one  can  see,  the  presence  of  the  hole 
has  reduced  the  peak  pressure  by  almost  6  GPa,  while  retarding  the  advance  of 
the  flame  front  by  2  cm.  Indeed,  even  at  t  =  72  ysec,  the  peak  pressure  is 
still  about  3  GPa  less  than  the  baseline  case  at  t  =  68  ysec. 

Figures  6  and  7  compare  the  effect  of  varying  the  width  (CQ)  of  the  crack 
while  keeping  the  length  constant  at  4  mm.  Figure  6  shows  the  pressure 
profiles  at  x  =  5.0  cm  for  various  values  of  CQ.  It  is  apparent  that  as  CQ 
increases,  the  profiles  converge.  Figure  7  shows  the  effect  of  the  various 
CQ's  on  the  flame  front  velocity.  As  can  be  seen  from  the  plots,  the  flame 
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front  achieves  a  constant  velocity  near  the  end  of  the  burn.  This  final 
Chapman-Jouget  velocity  is  one  of  the  critical  parameters  in  determining 
detonation.  As  one  can  see  from  the  plots,  the  final  velocity  decreases  as  CQ 

increases.  We  also  see  that  the  flame  profiles  again  converge  for  large 

values  of  C0. 

Figures  8  and  9  show  the  effect  of  varying  the  length  of  the  crack  while 
keeping  CQ  constant  at  a  value  of  0.1.  Figure  8  shows  the  pressure  versus 
distance  profiles  at  t  =  S8  usee  for  cracks  ranging  from  4  mm  to  16  mm  in 

length.  As  one  can  see,  varying  the  length  of  the  crack  has  a  more  severe 

effect  on  the  profile  than  varying  CQ.  For  the  case  where  W  -  0.16  (16  mm 
crack),  the  peak  pressure  in  the  bed  is  a  mere  fraction  of  the  value  obtained 
for  the  totally  confined  case.  Moreover,  Figure  9  shows  that  the  final  flame 
front  velocity  has  been  reduced  from  5.2  mm/usec  to  only  2.7  mm/ysec. 

While  results  such  as  Figure  5  shows  that  partial  confinement  does  indeed 
have  a  significant  effect  on  the  forming  detonation  wave,  one  is  still 
concerned  whether,  after  encountering  a  hole  in  the  confining  wall,  the 
reaction  wave  could  possibly  recover  and  still  form  a  detonation  wave  if  there 
were  a  sufficient  length  of  confined  bed  downstream  of  the  hole.  Therefore,  a 
condition  was  simulated  for  a  bed  of  20  cm  length  with  no  hole  in  the 
casing.  This  result  was  then  compared  to  another  20  cm  long  bed  in  which  a 
hole  with  W  =  0.02  (2  mm  in  length)  and  a  CQ  =  0.20  had  been  introduced  into 
the  casing.  Figure  10  shows  the  pressure  versus  distance  profile  for  the 
partially  confined  bed.  Also  shown  is  the  maximum  pressure  achieved  by  the 
totally  confined  bed.  As  one  can  see,  even  though  the  process  in  the 
partially  confined  case  required  10  usees  longer,  the  peak  pressure  in  the 
damaged  case  has  almost  completely  recovered  the  depressurization  suffered  due 
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to  the  hole.  Figure  11  further  confirms  this  recovery  by  showing  that  the 
final  velocities  of  the  two  flame  fronts  are  almost  equal  for  the  two  cases. 
This  result  was  somewhat  unexpected.  It  was  at  first  believed  that  the 
presence  of  such  a  hole  would  permanently  quench  a  detonation.  Of  course,  a 
crack  of  much  longer  width  would  likely  prevent  a  detonation  altogether. 

Finally,  a  case  was  run  in  which  the  conditions  for  the  crack  formation 
were  changed.  Instead  of  assuming  that  the  crack  already  existed  at 
prescribed  locations  at  t  *  0,  it  was  assumed  that  the  impact  with  the 
hardened  surface  had  instead  weakened  the  casing  so  that  if  the  pressure 
inside  the  bed  exceeded  a  prescribed  critical  pressure,  the  casing  would 
crack.  The  model  was  set  up  so  that  when  the  pressure  at  any  point  in  the  bed 
exceeded  5  kbar  (72,500  psia)  a  crack  with  CQ  =  0.1  and  W  =  0.01  would  form  at 
that  location. 

Figure  12  shows  the  pressure  versus  distance  profile  for  this  case.  As 
can  clearly  be  seen  the  pressure  in  the  damaged  explosive  never  exceeds  the 
prescribed  critical  pressure.  Therefore  it  would  be  impossible  for  a 
detonation  to  occur  if  such  physics  as  modeled  were  actually  valid.  Figure  13 
shows  the  flame  front  locus  for  the  totally  confined  case  and  for  two  cases 
where  the  casing  fails  at  a  prescribed  pressure  --  one  at  P  =  0.5  GPa  and  one 
at  P  =  1.0  GPa.  One  can  see  the  severe  reduction  in  the  final  flame  velocity 
caused  by  the  erupting  cracks. 


7.  CONCLUSIONS  AND  REMARKS 

The  results  presented  in  the  previous  section  verify  the  fact  that 
incomplete  confinement  of  the  explosive  media  does  have  a  strong  affect  on  the 


formation  of  a  steady-state  detonation  wave  in  a  fragmented  explosive.  The 
presence  of  a  hole  can  bring  about  a  very  significant  reduction  in  the  gas 
pressure  as  well  as  cause  a  severe  reduction  in  the  reaction  front  velocity. 
The  results  also  showed  that  if  a  sufficient  length  of  confined  bed  remained 
downstream  of  the  hole,  a  detonation  wave  could  recover  and  still  reach 
steady-state  detonation. 

The  procedure  outlined  in  this  section,  while  providing  results  which  are 
heuristical ly  reasonable,  does  make  some  major  simplifying  assumptions  in 
order  to  allow  us  to  use  the  existing,  proven  one-dimensional  flame  spreading 
model.  Perhaps  the  greatest  weakness  of  this  model  is  that  it  does  not 
simulate  the  actual  two  and/or  three  dimensional  processes  of  the  gas  and 
particles  flowing  out  of  the  hole.  Since  the  one-dimensional  model  allows 
only  flow  in  the  direction  of  the  deflagration  front,  it  is  impossible  to 
model  the  turning  of  the  gas  as  it  gains  y  and  z  velocity  components. 
Furthermore,  the  one-dimensional  model  does  not  allow  us  to  take  the  geometry 
of  the  bomb  or  the  cracks  into  account. 

Therefore,  we  thought  it  necessary  to  formulate  an  unsteady  flartie 
spreading  model  in  three  space  dimensions  to  allow  us  to  simulate  this  problem 
with  greater  accuracy.  In  such  a  three-dimensional  model,  the  hole  in  the 
confining  wall  would  simply  become  a  boundary  condition  which  allows  the 
normal  velocity  component  at  the  wall  to  be  nonzero.  (Rather  than  the  usual 
boundary  condition  which  states  that  the  normal  velocity  component  at  the  wall 
must  be  zero.) 
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Schematic  of  Impact  Damaged  Warhead 
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Figure  2.  Schematic  of  Explosive  Bed  Showing  Assumed  Crack  in  Casing 
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Figure  4.  Pressure  History  for  Baseline  Case  (no  crack  in  cas 
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lame  Front  Location  vs.  Time  for 
Various  Hole  Sizes 
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Figure  7.  Locus  of  Flame  Front  for  Various  Crack  Values  of  C 
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Flame  Front  Locus  for  Erupting  Cracks  for  Two  Different  Critical 
Pressures  (locus  of  cracks  is  also  shown) 


SECTION  III 


THE  THREE  DIMENSIONAL  MODEL 


1.  INTRODUCTION 

In  References  1-8,  the  separated  flow  approach  was  used  to  derive  the 
conservation  equations  for  an  unsteady,  two-phase  reactive  system  in  one  space 
dimension.  In  this  section,  the  same  type  of  analysis  is  followed  to  derive 
the  conservation  equations  for  a  three-dimensional  domain  in  order  to  properly 
model  the  dynamic  flow  processes  in  partially  contained  regions.  Each  phase 
will  be  considered  as  a  separate  fluid  flowing  through  its  own  control  volume, 
the  sum  of  the  two  volumes  being  equal  to  the  average  mixture  volume.  The 
conservation  equations  are  then  expressed  in  cylindrical  coordinates  with 
interphase  mass,  momentum  and  energy  transfer  terms  being  included  in  all 
three  independent  directions. 

Before  work  was  started  on  this  phase  of  the  project,  a  thorough  search 
of  the  fluid-mechanics  literature  was  made.  While  the  search  was  quite 
exhaustive,  very  little  work  on  transient  flows  was  found  which  could  be  of 
help  to  the  analysis  of  this  problem.  However,  a  paper  by  Markatos  and 
Kirkcaldy  (Reference  13),  was  found  in  which  a  numerical  model  based  on  a 
separated  flow  analysis  was  used  to  investigate  the  transient,  three- 
dimensional  reactive  flow  through  the  totally  confined  granulated  solid 
propellant  inside  a  gun  cartridge.  While  this  paper  did  not  present 
sufficient  details,  it  did  provide  a  useful  example  which  allowed  us  to 
compare  the  unsteady,  three-dimensional  conservation  equations  derived  in  this 
work  with  those  derived  independently  by  other  authors. 


2.  ASSUMPTIONS 


All  of  the  assumptions  stated  in  Section  II,  paragraph  2,  also  apply  for 
the  three-dimensional  case,  with  the  obvious  exceptions  of  assumptions  (b)  and  „ 
(e),  which  refer  to  the  one-dimensional  nature  of  the  equations.  It  will 
still  be  assumed  that  the  gas  escaping  from  the  domain  through  the  crack  in 
the  wall  will  leave  at.  the  choked  velocity.  However,  the  pseudo  one¬ 
dimensional  sink  term  (e)  used  to  approximate  the  lost  mass,  momentum,  and 
energy  will  not  apply  here.  Instead  the  choked  assumption  will  comprise  the 
boundary  condition  at  the  location  of  the  crack,  and  the  amount  of  mass, 
momentum,  and  energy  lost  will  be  found  directly  from  the  solution  of  the 
conservation  equations. 


3.  THE  UNSTEADY,  THREE-DIMENSIONAL  CONSERVATION  EQUATIONS 

The  unsteady  conservation  equations  are  presented  below  in  their  complete 
form.  A  detailed  derivation  of  these  equations,  based  on  the  separated  flow 
concept  is  presented  in  Appendix  C.  The  gas  and  particle  phase  momentum 
equations  in  the  r  ,  o  and  z  -  direction  are  presented  in  generic  form  for 
brevity,  with  Table  3  giving  the  specific  terms  for  each  phase  and  direction. 
The  conservation  equations  are: 
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Gas  Phase  Momenta  (3-equations) 
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Solid  Phase  Momenta  (3- equations) 
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where  *j  stands  for  ug,  vg,  cr  wg,  respectively  and  $2  stands  for  up,  vp,  or 

w„.  The  terms  S.  and  S.  represent  the  gas  and  solid  phase  source  terms, 

P  $1  *2 

respectively. 
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TABLE  3.  SOURCE  TERMS  FOR  MOMEMTUM  EQUATIONS 
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In  Equations  (18)-(20) 
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and  in  Equations  (21) - (23) 
P2  =  Pp  (!-♦) 


Gas  Phase  Energy 
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Solid  Phase  Energy 
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In  order  to  completely  describe  our  system  at  every  instant  of  time,  we 
need  to  know  the  instantaneous  values  of  og,  pp,  ug,  upt  vg,  vp,  wg,  wp,  Tg, 
Tpi  Pg,  Pp1  and  $  -  -  -  thirteen  variables.  Equations  ( 16) - ( 25 )  allow  us  to 
solve  for  ten  of  these  variables.  We  therefore  need  three  additional 
relations  to  achieve  closure.  These  equations  are: 
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(i)  an  equation  of  state  for  the  gas  phase 

(ii)  an  equation  of  state  for  the  solid  phase 

(iii)  a  material  stress-porosity  relation 

The  same  relations  used  in  the  analysis  described  in  Section  II  (and 
further  described  in  Appendix  A)  will  be  used  here  as  well.  Furthermore,  the 
same  relations  for  the  gas  mass  generation  r  and  the  heat  transfer  rate  Q,  as 
used  in  Section  II  will  be  used  here.  The  definition  of  the  drag  viscous 
interaction  will  be  assumed  to  hold  true  in  each  direction  independently, 
i  .e. , 
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In  comparing  Equations  ( 16) - (25)  to  the  one-dimensional  conservation 
equations  (Equations  (l)-(6)  in  Section  II),  one  can  see  that  the  main 
differences  (other  than  the  obvious  additional  terms  that  account  for  the 
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multi-dimensional  nature  of  the  new  equations)  are  the  addition  of  the 

term  in  the  r-direction  momentum  equations  and  the  — °  term  in  the 
©-direction  momentum  equations.  These  terms  are  the  radial  acceleration  and 
the  Coriolis  acceleration  terms,  respectively,  and  come  from  the  fact  that 
*  0  and  *  0  in  a  curvilinear  coordinate  system.  [Here,  r  is  the  unit 

do  do 

vector  in  the  radial  direction  and  e  is  the  unit  vector  in  the  azimuthal 
direction. ] 

Also  on  the  inspection  of  Equations  ( 18)- (23)  we  note  that  the  pressure 
gradient  terms  are  now  written  as  the  gradient  of  partial  pressure  (i.e. 

)  rather  than  the  gradient  of  pressure  (i.e.  ^  ),  as  is  sometimes 

found  in  the  literature  (Reference  13).  The  specific  form  of  this  term  has 
been  a  subject  of  debate  among  fluid  dynamicists  working  in  multi-phase 
flows.  We  feel  that  the  correct  form  should  be  ,  since  the  term 

represents  the  pressure  force  acting  on  the  control  volume  faces.  Since  we 
are  considering  each  phase  separately,  we  should  include  only  the  pressure 
acting  on  that  area  of  the  control  volume  face  which  is  occupied  by  the  gas 
phase.  Hence,  the  porosity  «  (gas  volume/total  volume)  must  be  included 
inside  the  derivative.  Naturally,  the  same  reasoning  holds  true  for  the  solid 
phase.  Therefore,  the  term  (l-$)  must  be  included  inside  the  derivative  of 
the  solid  phase  pressure  gradient  term. 

4.  DOMAIN 

The  first  approximation  to  the  actual  geometry  of  current  warhead 
containers  would  be  to  model  the  bomb  as  a  right  circular  cylinder,  as  is 
shown  in  Figure  14.  The  figure  shows  the  assumed  crack  in  the  outer  casing  as 
well  as  the  fragmented  explosive  inside.  To  be  able  to  model  realistic  crack 


configurations  one  would  have  to  include  the  entire  360°  cross  section  of  the 
cylinder  in  the  domain  of  the  model.  However,  one  of  the  difficulties  in 
tackling  this  problem  is  the  large  amount  of  computer  time  and  memory  which 
would  be  required  to  simulate  such  cracking.  Therefore,  in  an  effort  to 
reduce  the  size  of  the  domain,  it  was  decided  to  assume  a  certain  amount  of 
symmetry  in  setting  up  this  problem.  This  can  be  done  by  assuming  that  the 
damage  to  the  casing  occurs  in  such  a  manner  that  holes  form  in  a  symmetric 
ring  around  the  casing,  as  is  shown  in  Figure  15.  This  subdivides  the  domain 
of  the  cylindrical  bomb  into  symmetric  wedges,  the  flow  in  each  wedge  being 
identical  to  the  flow  in  all  the  other  wedges.  The  flow  inside  the  bomb  may 
then  be  completely  described  by  solving  the  flow  inside  a  single  wedge. 

Although  such  a  domain  limits  the  realism  of  this  problem  from  a  geometry 
standpoint,  this  assumption  still  retains  the  necessary  multi-dimensional 
features  of  the  transient  fluid  analysis.  Also,  it  is  not  currently  within 
the  scope  of  this  work  to  analyze  the  dynamics  of  the  casing  fracture,  but 
rather  to  develop  a  model  of  the  transient,  two-phase  reactive  fluid  mechanics 
involved  in  the  processes  of  detonation  transition  or  failure.  Clearly  the 
model  developed  herein,  under  this  specific  symmetry  assumption,  can  be 
extended  by  simply  adjusting  the  boundary  conditions  in  the  5  -direction  to 
include  the  full  domain,  once  a  super-computer  with  sufficient  speed  is 
avai lable. 

5.  BOUNDARY  CONDITIONS  AND  INITIAL  CONDITIONS 

Once  the  domain  of  the  model  has  been  fixed,  boundary  conditions  must  be 
properly  specified.  The  symmetric  domain  shown  in  Figure  15  has  six 
boundaries:  solid  walls  at  r  =  R,  z  -  0,  and  z  -  L,  and  symmetry  boundaries 
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at  r  -  0,  9  =  0|  and  e  =  02-  Obviously,  at  each  solid  wall  the  normal 
velocity  component  must  be  zero.  In  addition,  symmetry  allows  us  to  say  that 
the  radial  velocity  components  Ug  and  Up  must  be  zero  at  r  =  0,  and  all 
derivatives  with  respect  to  r  must  also  be  zero  here.  Symmetry  also  implies 
that  Vg  and  vp  must  be  zero  at  and  02  and  that  all  derivatives  with  respect 
to  9  must  be  zero  at  these  two  boundaries. 

As  mentioned  in  the  introduction  to  this  chapter,  the  velocity  at  the 
location  of  the  crack  will  be  assumed  to  be  choked.  This  assumption  forms  the 
last  boundary  condition  needed  to  solve  the  finite  difference  equations.  How 
these  boundary  conditions  are  incorporated  into  the  finite  difference  model  is 
discussed  in  detail  in  Section  IV,  paragraph  4. 

The  bed  of  explosive  particles  is  assumed  to  be  initially  quiescent  and 
at  a  prescribed  nominal  temperature  and  pressure.  Combustion  is  initiated  at 
time  t  =  0  by  assuming  that  a  hot  burning  zone  exists  at  one  end  of  the 
domain,  with  the  temperature  in  this  zone  enough  above  the  explosive  ignition 
temperature  to  ignite  the  particles.  The  initial  value  input  is  given  in 
Section  V. 


3-DIMENSIONAL 
WEDGE  MODEL 


SECTION  IV 


THE  FINITE  DIFFERENCE  EQUATIONS 


1.  INTRODUCTION 

The  unsteady  two-phase  conservation  equations  presented  in  Section  III 
form  a  set  of  coupled,  non-linear,  hyperbolic  partial  differential 
equations.  As  was  the  case  in  the  pseudo  two-dimensional  formulation 
presented  in  Section  II,  these  equations  also  will  not  have  analytical 
solutions.  Futhermore,  a  Method  of  Lines  (MOL)  approach  such  as  the  one  used 
in  solving  the  one-dimensional  equations  would  not  be  feasible  for  the  three- 
dimensional  case  due  to  the  great  increase  in  the  size  of  the  domain.  The 
high-order  ODE  solvers  used  in  the  Method  of  Lines  require  a  large  number  of 
calculations  to  be  performed  when  solving  the  unsteady  terms  in  the 
conservation  equations.  Therefore,  while  such  an  approach  may  be  feasible  for 
the  domain  of  the  one-dimensional  case,  it  would  prove  to  be  cost  prohibitive 
for  the  much  larger  three-dimensional  domain.  Instead,  a  new  integration 
method  was  developed  as  part  of  this  research. 

In  this  section,  we  will  detail  the  development  of  a  new  numerical  scheme 
which  will  be  used  to  integrate  Equations  (16)-(25).  The  discretization  of 
the  domain  into  a  mesh  of  grid  points  is  discussed  in  paragraph  2,  while  the 
actual  finite  difference  scheme  is  discussed  in  paragraph  3.  Paragraph  4 
details  the  finite  differencing  of  the  boundary  conditions.  Finally, 
stability  considerations  for  the  finite  difference  scheme  are  discussed  in 
paragraph  5. 


2.  THE  FINITE  DIFFERENCE  GRID 


In  applying  a  finite  difference  method  to  the  solution  of  any  system  of 
equations,  the  continuous  domain  of  the  problem  must  be  replaced  by  a  mesh  of 
grid  nodes,  each  node  representing  a  location  where  the  values  of  certain 
dependent  variables  are  stored.  In  this  analysis,  a  staggered  grid  was  used 
to  represent  the  domain.  Figure  16  shows  a  single  cylindrical,  three- 
dimensional  control  volume  element  employing  the  staggered  grid.  All  scalar 
variables  are  stored  in  the  node  at  the  center  of  the  control  volume,  while 
the  bulk  momentums  in  the  r  -direction  (  p.u  and  p^u  )  are  stored  at  each  of 

'  i  g  c  p' 

the  nodes  in  the  r  -faces  of  the  control  volume,  the  bulk  momentum  in  the  e  - 

direction  (  p.v  and  p,v  1  are  stored  at  each  of  the  nodes  in  the  9  -faces, 

1  9  ^  P 

and  the  bulk  momentums  in  the  z  -direction  (  p^  and  p2wp)  are  stored  at 
each  of  the  nodes  in  the  z  -faces. 

There  are  several  reasons  why  it  is  advantageous  to  use  the  staggered 
grid  instead  of  the  conventional  grid.  First,  as  will  be  shown  in  paragraph 
4,  it  makes  the  finite  differencing  of  the  boundary  conditions  considerably 
easier.  It  also  enables  one  to  get  second-order  accurate  centered  space 
differences  while  having  to  span  only  half  the  distance  required  for  taking 
such  a  difference  on  a  conventional  grid.  The  advantage  of  this  can  be 
illustrated  by  the  following  one-dimensional  example.  Consider  the  two  grids 
shown  in  the  Figure  17.  The  upper  grid  is  a  conventional  grid,  while  the 
staggered  grid  is  shown  below  that. 
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Now  assume  one  wishes  to  evaluate  the  pressure  gradient  term  in  the 
momentum  equation.  On  the  conventional  grid,  this  term  is  evaluated  by  second 
order  centered  difference  as: 


aJP 

3X 


Pi_1)/2AX 


If  we  apply  this  equation  at  node  i=2,  we  get  for  our  example, 
£  =  (100-100)/2ax  =  0 


This  obviously  cannot  be  true  given  the  pressure  field  depicted  in  the  figure. 

If  we  now  evaluate  this  finite  difference  on  the  staggered  grid,  we  only 
need  to  take  the  difference  over  a  total  step  of  ax  (rather  than  2ax  as 
before),  since  values  of  Pg  are  available  at  nodes  that  are  a  distance  of 
ax/2  and  -  ax/2  from  the  velocity  node  (where  the  momentum  equation  will  be 
evaluated).  Thus  at  i-2. 


^  -  <**i*  %  -  <■,_  *>/“ 


or 


3Pn 

>  (100-50) /ax  -  50/ax 

o  X 


As  one  can  see,  the  scheme  using  the  staggered  grid  yields  a  much  more 
accurate  result  for  this  typical  example. 
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Figure  18  shows  the  r,  9  -plane  and  r,  z-plane  projections  of  our  wedge- 
shaped  domain  with  the  staggered  grid  superimposed  over  it.  The  continuity 
equation  (Equations  (16)  and  (17)),  which  are  used  to  solve  for  Pl  and  p2  , 
and  the  energy  equations  (Equations  (24)  and  (25)),  which  are  used  to  solve 
for  Egt  and  Ept,  will  be  solved  at  the  scalar  nodes  (denoted  by  the  symbol  + 
in  Figure  18)  since  these  are  scalar  quantities.  The  r  -momentum  equations 
(Equations  (18)  and  (19))  will  be  solved  at  the  r  -momentum  nodes  (denoted  by 
the  symbol  ),  the  I  -momentum  equations  (Equations  (20)  and  (21))  will  be 
solved  at  the  $  -momentum  nodes  (denoted  by  the  symbol  ),  and  the  z  - 
momentum  equations  (Equations  (22)  and  (23))  will  be  solved  at  the  z  -momentum 
nodes  (denoted  by  the  symbol  ). 

When  employing  the  staggered  grid,  the  situation  may  arise  where  the 
value  of  a  variable  is  needed  at  a  location  where  it  is  not  defined.  In  this 
case  a  simple  average  will  be  taken  to  find  the  value  at  the  desired  point. 

For  example,  suppose  one  needs  to  know  the  value  of  the  velocity,  u,  at  the 
scalar  node  i=l  shown  in  the  one-dimensional  depiction  (Figure  19). 

Since  the  u  values  are  not  stored  at  i=l,  we  must  average  the  values  stored  at 
i=  %  and  i=l  *s,  together,  to  form  u._^.  Thus 

“i  ■  is  5/2 

Although  the  necessity  for  such  averaging  may  complicate  the  numerical 
scheme  somewhat,  it  is  actually  advantageous  as  far  as  the  accuracy  of  the 
calculations  is  concerned.  The  averaged  value  is  influenced  by  the  value  of 
the  variable  at  two  nodes,  while  the  variable  which  is  read  directly  from  a 
location  where  its  value  is  known  is  influenced  only  by  a  single  value.  In 
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other  words,  the  averaging  operation  passes  information  about  gradients  as 
well  as  providing  the  value  needed. 

Finally,  one  will  note  that  in  our  scheme,  momentum  rather  than  velocity 
is  stored  at  the  nodes  located  in  the  control  volume  faces.  This  was  done 
mainly  to  reduce  the  amount  of  computer  storage  required  for  the  code. 

3.  THE  FINITE  DIFFERENCE  EQUATIONS 

Once  the  domain  has  been  discretized,  we  must  next  cast  the  continuous 
conservation  equations  (16)-(25)  into  finite  difference  form.  We  will 
basically  be  using  the  leapfrog  scheme  Reference  14,  which  means  that  the 
unsteady  term  will  be  evaluated  at  time  levels  t  +  at  and  t  -  At  ,  while  the 
flux  and  source/sink  terms  on  the  right  hand  sides  of  equations  (16)-(25)  will 
be  evaluated  at  the  current  time,  t. 

The  leapfrog  scheme  was  chosen  after  we  reviewed  a  paper  by  Williams 
(Reference  15)  in  which  he  successfully  used  the  leapfrog  scheme  to  integrate 
the  unsteady  three-dimensional  conservation  equations  for  an  incompressible, 
single-phase  viscous  fluid.  The  similarity  between  the  equations  integrated 
in  Williams  paper  and  the  equations  handled  in  this  work  prompted  us  to 
consider  a  leapfrog-type  scheme.  While  the  compressible,  two-phase  reactive 
nature  of  the  flow  dealt  with  in  this  paper  resulted  in  some  major  differences 
between  the  numerical  scheme  developed  by  Williams  and  the  scheme  developed 
herein,  the  usefulness  of  the  work  done  in  Reference  15  in  aiding  us  in  our 
analysis  is  still  acknowledged. 

Consider  the  following  shorthand  notation  (where  x  is  a  generalized 
independent  variable  and  q  is  a  generalized  dependent  variable): 
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sxq  =  [q  (x+ax/2)  -  q  (x-ax/2)]/ax 


(30) 


q*  5  [q  (x+ax/2)  +  q  (x-ax/2)]/2 


(31) 


Equation  (30)  is  a  difference  operator  while  Equation  (31)  defines  an 
averaging  operator.  Furthermore,  let  us  define 

m  s  o . u  -  qas  bulk  momentum  in  the  r-direction 

r  in  J 


m  e  p2u  -  particle  bulk  momentum  in  the  r-direction 
r  2  P 


m  =  oiV  -  gas  bulk  momentum  in  the  e-direction 
Si  1  g 


m  =  0_v  -  particle  bulk  momentum  in  the  ^-direction 
e2  2  p 


m  =  p,w  -  gas  bulk  momentum  in  the  z-direction 
2 1  9 


m  =  o2w  -  particle  bulk  momentum  in  the  z-direction 
^  2  P 


We  can  then  write  Equations  (16)-(25)  based  on  our  staggered  grid  in 
finite  difference  form  as: 

Gas  Continuity  (applied  at  scalar  nodes) 


.  -t 

v> = 


•F‘A,r>  -  W  *  r 


(32) 
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Particle  Continuity  (applied  at  scalar  nodes) 


-  -  7  MV'  -  7  se<me;>  '  M"*,1  '  r  <33> 

6as  r  -Direction  Momentum  (applied  at  r  -momentum  nodes) 

5tS5,=  -  r  4rI''(^1)Vo1]  -  7  S9  [(ij,  (",,/?!)  I 

-  6  [m~z(m  h\)  +  (m  )2/G\r) 

*•  1  i  4 1  1 1 

-  UP,)  -  dV  /p^  -  m  h\) 

r  r  1  r2 

+  r  m  hr2  (34) 

1  2 

Particle  r  -Direction  Momentum  (applied  at  r  -momentum  nodes) 

Vr,3 ' 7  Vtr<5r//o'd  -  7  6eK!(”e/»°)l 

"6z("'r2(m22/^)l  + 

-  «r(P2)  +  D  (mr  h\-  m rhr2) 

-  Tr  m  fpr2  (35) 

r  2 

Gas  e  -Direction  Momentum  (applied  at  §  -momentum  node) 
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SA-  -  7  srK/"V/‘“‘)  I  -  7 


K<P.>  -  D  <V»!  -  ma/^> 


+  79 


Particle  $  -Direction  Momentum  (applied  at  5  -momentum  nodes) 


Ve,  -  -  7  VK  -  7  «e[(»n)‘/»2l 


1  *  r/cr-S\2 


-  {2l’n82(",2/"1>l  -  W°2r) 

-  7  «8(Pi)  +  D  (ly/^  -  i"9j/o|) 


— e  e 

r  »e/»2 


fias  I  -Direction  Momentum  (applied  at  z  -momentum  nodes) 


8tmz -  7  VK 


r[^,(”r/^  i  'i 


49 


(applied  at  scalar  nodes) 


Particle  Energy 


6tE?  =  -  7  vK  r+  (p>I)rH 

9+ 

-  Sz{m22  [E77p2  +  (P>7)Z|} 

+Q+Dm^/p2  +  0  i*  /p2 
”2  y  2 

+  oiz/».  +  r  l^hem  -  (^/.^Vz 

-  (m.  /p2)  /2  -  (m  /p2)  /2]  (41) 

0  2  *2 

In  these  equations. 


0  =  (ug/4rpZ)  • 


4.  THE  BOUNDARY  CONDITIONS 


Once  the  domain  has  been  discretized  and  an  integration  scheme  has  been 


devised,  we  must  next  turn  to  implementing  the  boundary  conditions.  As 
mentioned  earlier,  there  are  six  different  boundaries  where  tne  equations  must 
be  satisfied.  We  will  consider  each  boundary  separately. 


(1)  r  =  R 


Here  we  have  a  solid  wall,  which  means  that  the  normal  velocity 
components  (Ug  and  Up)  must  be  zero.  Figure  18  shows  that  at  this  boundary, 
we  have  only  r  -momentum  nodes.  Since  only  the  ?  -momentum  equations  are 
evaluated  at  these  nodes,  this  boundary  condition  can  easily  be  met  by 
initializing  the  value  of  the  r  -momentum  variables  stored  at  these  nodes  to 
zero  and  ordering  the  code  to  skip  calculating  the  change  in  the  r  -momentum 
(i.e.  skip  evaluating  Equations  (34)  and  (35))  when  r  =  R.  This  will  cause 
the  values  of  the  r  -momentum  to  remain  at  zero,  thus  satisfying  the  boundary 
condition. 

One  must  also  be  able  to  evaluate  Equations  (32),  (33)  and  (36)-(41)  at 
nodes  which  are  not  directly  on  the  boundary,  but  are  a  distance  of  &r/2  in 
from  the  wall  at  r  *  R.  The  evaluation  of  these  equations  at  these  particular 
nodes  is  affected  by  the  wall  boundary  in  the  evaluation  of  the  r  -direction 
flu*  term  (i.e.  the  first  term  on  the  right  hand  side  in  Equations  ( 32)— 

(41)). 

The  problem  and  its  solution  are  perhaps  best  described  by  an  example. 
Consider  attempting  to  evaluate  the  r  -Hux  term  in  Equation  (40)  (i.e. 

5  (m  [  ( E t /d +  ( p i /p i T ] } )  at  the  scalar  node  adjacent  to  the  boundary  at 

r  r  i 

r  =  R  as  shown  in  Figure  20.  Evaluating  the  flux  term  as  described  by 
Equation  (30),  we  can  see  that  we  must  find  the  values  of  E^c,  and  P1/p1  at 
i  +  5j  and  i  -  Jj.  However,  since  E,,  o,,  and  are  all  scalar  quantities, 
their  values  are  not  known  at  these  points  and  must  therefore  be  found  by 
interpolation  as  the  averaging  operators  in  the  flux  term  indicate.  This 
creates  an  apparent  problem  when  one  tries  to  average  the  values  at  (i)  and 
(i+1)  to  get  the  value  at  < i^4*) »  since  (i+1)  lies  outside  the  domain. 
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However,  one  will  note  that  the  quantity  (Et/p,  +  pj/pJ  is  always 
multiplied  by  the  momentum  in  the  r  -direction  to  find  the  flux,  and  since 
the  f  -momentum  at  the  wall  at  r  =  R  is  always  zero,  these  variables  can  have 
any  finite  values  at  (i+1).  Once  the  average  has  been  taken  and  the  result 
multiplied  by  the  ?  -momentum  at  (1+**) (which  is,  of  course,  zero),  the  flux 
through  the  wall  in  the  ?  -direction  will  have  the  correct  value  of  zero. 

This  is  why  these  nodes  are  refered  to  as  false  exterior  nodes.  They  lie 
outside  the  domain  and  are  initialized  to  some  arbitrary  value,  at  which  they 
remain  throughout  the  run.  They  are  used  purely  to  fulfill  the  requirements 
of  the  numerical  method  and  the  griding  scheme. 

(2)  r  -  0 

The  symmetry  boundary  conditions  requires  that  up  and  ug  must  be  equal  to 
zero  at  r  *  0  and  that  all  derivatives  with  respect  to  r  must  be  equal  to 
zero.  Again,  we  can  note  from  Figure  21  shown  below  that  there  are  only  r  - 
momentum  nodes  at  r  =  0,  The  first  boundary  condition  is  again  met  by 
initializing  the  values  of  the  gas  and  particle  phase'momenta  at  these  nodes 
to  zero  and  instructing  the  code  to  skip  solving  the  r  -momentum  equations 
whenever  r  =  0.  Note  that  this  procedure  also  avoids  the  singularity  that 
occurs  in  the  1/r  terms  as  r+Q  .  Since  we  have  only  r  -momentum  nodes  at  r* 0, 
and  since  the  code  skips  solving  these  equations  at  that  point,  the 
singularity  never  arises. 

The  symmetry  condition  can  be  met  by  using  reflection  points  about  the  z- 
axis  at  r=0,  and  continuously  updating  the  values  of  the  variables  stored  at 
these  nodes  to  match  the  nodes  directly  across  the  z-axis.  However,  referring 
to  Figure  21,  we  can  see  that  again  the  only  interaction  between  the  equations 
evaluated  at  nodes  which  are  a  distance  of  Ar/Z  from  the  boundary  at  r=0  and 


>  u* 

1  ■  -  ■■  a .  a  ■  a  ■  \ 


J»  -  *  ■  - 
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the  symmetry  boundary  at  r=0  is  in  the  r  -direction  flux  term,  just  as  in 
paragraph  Performing  the  same  analysis  which  we  performed  in  paragraph 

(1),  we  again  see  that  the  values  of  the  variables  at  these  symmetry  nodes  can 
be  completely  arbitrary.  Once  the  averaging  is  performed  and  the  result  is 
multiplied  by  the  r  -momentum  (which  is,  of  course,  zero),  the  flux  through 
the  boundary  at  r=0  will  have  its  proper  value  of  zero.  It  is  therefore  not 
strictly  necessary  that  the  values  at  the  reflection  points  be  continuously 
updated  to  enforce  the  symmetry  condition,  although  it  would  be  perfectly 
proper  to  do  so.  Note  that  the  only  reason  that  this  is  possible  is  because 
we  are  using  a  staggered  grid  and  our  equations  have  no  space  derivatives 
higher  than  first  order. 

(3)  z=0  (solid  wall) 

Here  we  again  have  a  solid  wall,  leading  to  the  boundary  condition  that 
the  velocity  components  wg  and  wp  must  be  zero.  This  boundary  condition  is 
imposed  in  the  same  manner  that  the  boundary  conditions  on  the  radial  velocity 
components  Ug  and  up  were  imposed  in  paragraph  (1).  Again  false  exterior 
nodes  are  used  to  facilitate  the  numerical  scheme  at  the  nodes  which  are  at  a 
distance  of  &z/2  from  the  boundary.  Figure  22  shows  the  details. 

(4)  z  =  L  (far  wall) 

The  same  boundary  conditions  that  applied  in  (3)  apply  here 

(5)  8  =  Q  j  (wedge  symmetry) 

Here  we  again  utilize  symmetry  with  the  boundary  condition  being  that 
tangential  velocity  components,  Vg  and  vp  at  these  nodes  are  zero.  This 
condition  is  enforced  in  the  same  manner  that  the  symmetry  condition  in 
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paragraph  (2)  is  enforced.  (See  Figure  23  for  details.) 

(6)  9  =  e2  (wedge  symmetry) 

The  same  conditions  that  applied  in  paragraph  (5)  apply  here. 

Open  vent  (r  =  R) 

Finally,  we  must  consider  the  boundary  conditions  at  r  =  R  at  the 
location  of  the  prescribed  holes  in  the  outer  wall.  Since  the  fluid  and 
particle  momentum  at  the  wall  in  this  case  will  no  longer  be  zero,  the  scheme 
which  was  devised  in  paragraph  (1)  will  not  apply.  Instead,  we  will  have  to 
use  one-sided  interpolations  to  determine  the  necessary  values  at  r  =  R. 

As  stated  before,  the  gas  will  be  assumed  to  be  leaving  the  hole  at  the 
choked  gas  velocity,  which  is  a  function  of  the  local  pressure  and  temperature 
below  the  hole.  Choked  flow  means  that  the  gas  velocity  is  equal  to  the  local 
sound  speed.  The  speed  of  sound  in  the  non-ideal  gas  we  are  using  here,  as  is 
derived  in  Appendix  B,  is  given  by: 


a  = 


C  (l+2blP  )  +  R  (l+blP!)  1 

q _ _ _ _  p  i  2 

Cv  pg(1+b>°g)  9 

9  s  y 


(8.9) 


In  addition  to  the  speed  of  sound,  it  will  be  necessary  to  know  the  value 
of  several  other  variables  of  the  fluid  at  the  wall.  These  may  be  found  via  a 
Taylor  Series  expansion  about  the  node  nearest  the  wall.  Referring  to  Figure 
24,  the  general  dependant  variable  W  defined  at  location  (i)  may  be 
approximated  at  (i+H)  as: 
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(42) 


W 


i+H 


w  +  AI  iW 
wi  2  ar 


Since  we  are  adjacent  to  the  wall,  the  derivative  cannot  be 

d  r 

approximated  as  a  centered  second-order  finite  difference.  Instead,  we  can 
write  this  term  as  a  one-sided  second  order  finite  difference: 


aW 

ar 


<r 


3W.-  4W._l+ 


i-2 


2hr 


(43) 


Therefore,  our  approximation  becomes 


W 


7V 


4W 


i-I 


W. 

V 


(44) 


If  the  hole  is  assumed  to  be  larger  than  one  node  in  the  ?-  or  §  -direction, 
it  is  assumed  that  the  velocity  at  the  node  at  the  center  of  the  hole  is 
choked,  and  that  the  velocities  at  the  remaining  hole  nodes  constitute  a 
linear  profile  so  that  the  u  velocity  is  zero  at  the  edge  of  the  hole.  (See 
Figure  25) 

The  velocity  of  the  particles  leaving  the  domain  through  the  hole  is 
evaluated  by  calculating  the  drag  force  on  the  particles  at  the  hole  location 


at  B  ’’0 

where  F,,  =  [)  (2*r  2) 
u  r  p 


(45) 


then 


3  (p2up)  55  F0dt 


(46) 


*>.  \  »,  *- 
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Integrating  both  sides,  we  get: 

/t+&td(p2un)  =  /t+atFndt  (47) 

t-At  P  t-At  U 

Since  we  assume  in  our  discretization  that  all  properties  remain  constant  over 
the  time  step  At,  FQ*  Fq( t)  .  Then 

/t+atd(p,u  )  =  Fn  ;t+Atdt  (48) 

t-At  ‘  P  U  t-At 
or 

(o2up)tMt-  (ojUp)1’41  -  f„.24t  (49) 

and  finally, 

<p,up)tMt  .  FD.  24t  *  (oaup)t-fit  (50) 

Once  a  value  for  the  particle  momentum  at  the  hole  is  known,  the  rest  of 
the  terms  needed  for  evaluating  the  numerical  scheme  can  be  found  via  the  one¬ 
sided  interpolation  described  above. 

5.  STABILITY 

The  numerical  solution  of  partial  differential  equations  by  finite 
difference  methods  always  presents  a  problem  called  stability.  A  finite 
difference  scheme,  while  providing  an  accurate  representation  of  the  original 
analytic  equation  and  boundary  conditions,  may  still  yield  unsatisfactory 
results  due  to  oscillations  and  explosive  growth  of  the  output  caused  by 
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instabilities  inherent  to  the  differencing  scheme  being  used.  Most  textbooks 
on  numerical  methods  treat  this  topic  in  great  detail.  In  this  section,  we 
will  look  at  several  sources  of  instability  inherent  to  the  scheme  presented 
in  paragraph  3  and  how  they  can  be  controlled. 

I 

In  any  explicit  differencing  scheme,  there  exists  a  bound  on  the  maximum 
size  of  the  time  step.  Our  conservation  equations  can  be  shown  to  be 
hyperbolic  partial  differential  equations.  If  the  size  of  the  time  step 

i 

exceeds  this  bound,  oscillations  will  set  in  and  the  output  will  soon  become 
unstable.  This  limiting  value  can  be  calculated  via  the  Von  Neumann  approach 

(see,  for  instance,  Richtmeyer  and  Morton,  Reference  16).  However,  for  the 

i 

system  of  ten  coupled,  nonlinear  equations  with  source  terms  which  form  our 

finite  difference  scheme,  the  evaluation  of  the  Von  Neumann  stability  criteria 

would  in  itself  constitute  a  formidable  task.  It  was  therefore  decided  to 

I 

defer  calculating  the  optimum  time  Step  and  simply  run  the  model  with  a 
suitable  constant  time  step,  at  ,  which  must  be  determined  by  trial  and 
error.  Once  the  suitability  of  the  scheme  has  been  proven,  then  a 
concentrated  effort  can  be  made  to  determine  the  formula  for  the  optimum 
allowable  time  step. 

A  second  source  of  instability  is  caused  by  the  fact  that  centered  time 
differencing,  rather  than  forward  time  differencing  was  used  to  represent  the 
unsteady  terms  in  Equations  (16)-(25).  Centered  time  differencing  was  chosen 
because  Kurihara  (Reference  14)  has  shown  that  it  provides  less  damping  of  the 
kinetic  energy  than  other  methods  of  time  differencing.  The  centered-time 
difference  scheme  can  be  represented  as: 

(Wt+1-  Wt-l)/2At  =  F  (r.e.z)1  +  s*  (51) 

$ 
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where  W  is  some  generalized  conserved  quantity,  F(r,o,z)  represents  the  flux 
terms  and  represents  the  source  terms.  Solving  for  we  get 

Wt+1  =  2at  (FCr.e.z)1  +  S  +  Wt_1  (52) 

Since  the  first  term  on  the  right-hand  side  is  added  to  the  value  of  W  at  time 
level  t-1,  there  is  a  possibility  that  W  at  time  level  t+1  may  be  less  than 
the  value  of  W  at  time  level  t,  even  if  2at  [F(r,e,z)t  +  S  t]  is  a  positive 
quantity.  Once  this  occurs,  the  solution  begins  to  oscillate  until  eventually 
separate  solutions  form  at  odd  and  even  time  steps. 

This  oscillation  is  controlled  by  adding  in  a  step  which  averages  the 
conserved  variables  over  adjacent  time  steps,  as  was  suggested  by  Williams 
(Reference  15).  Thus,  one  writes 

W*t_l  =  (W1  +  WtA)/Z  (53) 

r 

In  this  work,  the  variables  were  averaged  at  each  adjacent  time  step  to 
control  this  oscillation.  It  was  therefore  decided  that  a  time  step  of  3/2  at 
rather  than  ?at  should  be  used,  since  this  averaging  process  effectively 
reduces  the  size  of  the  time  step.  Thus,  we  assume 

Wt+1  =  (3/2)  at-F(r,e,z)t  +  W*1'1  (54) 


To  show  the  effect  of  this  procedure  it  will  be  useful  to  consider  one  of 
the  computational  results.  Figure  26  shows  a  comparison  of  the  flame  front 
locus  versus  time  using  time  steps  of  3/2at  and  2at,  respectively.  This  is 
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compared  with  the  flame  front  history  from  the  one-dimensional  model  described 
in  Section  II.  In  all  cases  the  walls  are  impermeable.  The  three-dimensional 
calculation  was  in  this  case  uniformally  initiated  to  make  the  comparisons 
meaningful.  As  can  be  seen,  the  use  of  3/2at  as  the  time  step  yields  results 
which  are  much  closer  to  the  results  from  the  previous  one-dimensional  model, 
a  result  which  is  fairly  accurate. 

Another  source  of  instability  results  from  the  nature  of  the  equations 
being  differentiated.  Equations  ( 16) - ( 25 )  represent  a  set  of  nonlinear, 
hyperbolic  equations.  The  numerical  integration  of  such  equations  can  become 
unstable  due  to  explosive  growth  of  the  total  energy  of  the  system.  This 
instability  is  caused  by  aliasing,  in  which  waves  that  are  too  short  to  be 
resolved  by  the  grid  are  misinterpreted  as  waves  of  longer  wavelength.  Figure 
27  illustrates  this  phenomenon. 

The  solid  line  in  Figure  27  denotes  the  actual  waveform.  However,  since 
this  continuous  wave  has  been  discretized,  the  code  senses  only  the  amplitude 
at  each  node.  It  therefore  incorrectly  interprets  the  wave  as  having  a 
waveform  shown  by  the  dashed  line.  Miyakoda  (Reference  17)  has  shown  that 
high  frequency  waves  produced  by  convective  terms  in  the  governing  equations 
can  in  this  manner  produce  non-physical  increases  in  the  total  energy  of  the 
system.  However,  Arakawa  (Reference  18)  has  shown  that  if  the  convective 
terms  in  the  total  derivative  form,  as  they  are  in  Equations  (32)-(41), 
aliasing  will  be  controlled.  The  preceeding  discussion  on  aliasing  and  its 
control  is  condensed  from  the  paper  by  Williams  (Reference  15)  and  the  reader 
is  referred  there  for  further  information  on  the  subject.  More  on  aliasing 
can  be  found  in  the  notes  by  Wilhelmson  (Reference  19). 

Another  common  source  of  instability  in  the  finite  differencing  of 
inviscid  hyperbolic  equations  stems  from  the  formation  of  shock  waves  in  the 
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computational  domain.  The  presence  of  a  shock  can  cause  near  infinite 
gradients  to  occur  at  the  front  of  the  wave.  These  steep  gradients  quickly 
cause  the  governing  equations  to  become  unstable. 

This  type  of  instability  is  caused  by  the  fact  that  the  second  derivative 
viscous  terms,  which  are  found  in  the  complete  Navier-Stokes  equations,  are 
not  included  in  the  pseudo-inviscid  equations  used  here.  The  viscous  nature 
of  the  real  flow  would  tend  to  smear  the  shock  discontinuity  over  a  finite 
(albeit  small)  distance,  thus  decreasing  the  severity  of  the  gradient  across 
the  front. 

Control  of  this  instability  in  an  inviscid  set  of  equations  can  be 
achieved  by  adding  viscous-like  artificial  diffusion  terms  to  the  right  hand 
side  of  Equations  (32)-(4l)  which  take  the  form  of  second  derivatives.  The 
text  by  Ames  (Reference  20)  presents  sufficient  background.  These  terms  act 
to  smooth  out  the  sharp  discontinuities  which  occur  when  a  shock  forms  in  the 
domain.  The  artificial  viscosity  term  used  in  this  work  is  similar  to  that 
used  in  Hyman's  Predictor-Corrector  Method  described  by  Sod  (Reference  21). 

r 

For  the  explicit  scheme  used  here,  it  was  necessary  to  lag  the  evaluation  of 
the  artificial  viscosity  terms  (i.e.  evaluate  them  based  on  the  values  of  the 
variables  at  the  previous  time  step)  to  ensure  stability  (Wilhelmson,  class 
notes  Reference  19).  Since  no  journal  articles  were  found  which  described  the 
application  of  artificial  viscosity  to  three-dimensional  problems,  one- 
dimensional  strategies  were  used  in  applying  these  diffusion  terms  to 
Equations  (32) - (41) .  One-dimensional  artificial  viscosity  terms  were  added  to 
the  equations  where  numerical  experiments  showed  them  to  be  necessary. 

During  the  development  of  the  computer  code,  it  was  found  that  the 
solutions  to  the  momentum  equations  quickly  became  unstable  unless  an 
unacceptably  small  time  step  was  used  (typically  10~®  seconds).  The  root  of 


this  instability  was  round  to  be  in  the  interaction  between  the  pressure 
gradient  and  drag  terms.  As  the  pressure  front  moves  through  undisturbed 
regions  of  the  bed,  the  strong  gradients  across  the  front  induce  rapid 
increases  in  the  gas  velocity.  The  difference  between  the  gas  and  particle 
velocities  would  then  give  rise  to  strong  drag  which  acts  to  retard  the  motion  - 
of  the  gas.  Realistically,  the  drag  should  continuously  adjust  itself  as  the 
relative  velocity  between  the  phases  changes.  However,  since  we  have 
discretized  this  process,  the  model  assumes  that  the  drag  remaining  constant 
over  the  entire  time  step.  If  the  drag  is  large  enough,  this  assumption  gives 
rise  to  the  situation  where  the  drag  not  only  causes  the  flow  to  slow  down, 
but  actually  reverses  its  direction.  This  is  physically  impossible,  since 
once  the  relative  velocity  between  the  phases  becomes  zero,  the  drag  would  be 
zero  and  there  would  be  no  force  to  retard  the  flow  any  further  and  cause  it 
to  become  negative. 

Once  the  gas  velocity  becomes  negative,  the  drag  begins  acting  in  the 
same  direction  as  the  pressure  gradient,  causing  the  velocity  to  again 
increase  sharply,  again  becoming  positive.  This  in  turn  again  causes  large 
drag  values  which  cause  the  gas  velocity  to  become  negative.  This  oscillation 
between  positive  and  negative  velocities  at  adjacent  time  steps  grows  until 
the  scheme  becomes  completely  unstable. 

This  problem  was  solved  by  using  a  predictor/corrector  strategy  when 
evaluating  the  momentum  equations  (Equations  (34)-(39)).  For  each  direction 
(radial,  azimenthal,  and  axial),  the  solid  and  gas  phase  momentum  equations 
were  solved  exactly  as  shown  in  Equations  ( 34 ) - ( 39 ) ,  with  all  of  the  terms  on 
the  right-hand  side  being  evaluated  at  the  current  time  level,  t.  This 
predicted  a  temporary  new  value  for  the  momenta  at  the  next  time  level,  t+1. 
These  temporary  values  were  then  averaged  with  the  current  values  of  the  gas 
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and  solid  phase  momentum  to  give  the  average  momentum  over  the  time  step. 

At.  Next  the  average  phase  velocities  over  the  time  step  were  found  by 
dividing  out  the  appropriate  phase  density.  Finally,  the  corrector  step  was 
performed  by  re-evaluating  the  momentum  equation  for  each  phase,  again 
evaluating  the  right-hand  side  at  the  current  time  level  (just  as  in  the 
predictor  step)  except  for  the  drag,  which  was  evaluated  based  on  the  average 


phase  velocities.  The  strategy  in  schematic  form  is,  then: 


(l)  OUiX,  =  F-<r’e'2)i,;i,k  "  +  <55> 


WTL  ■  F!(r,9,2)i,J,k  *  S*2  *  <56> 


(31  l(W,)^;kMW,)2J>k|/3  (57) 


*  tM*)^.j.k!/2  (58) 


(3)  (vi>;,j,k  ■  tf>;.Jfk  /  .f 


■  tf>;,jik/  .5 


(59) 


(60) 


<4>  '  Fi(r’e»z^i,j,k+  C-  °l^)*,J,k-  (^.J.kl 


(61) 


(W2) 


t+1 

i  J.k 


(62) 


Here  and  W£  represent  p^Ug  and  P?up»  Pivg  and  p£vp»  or  pjwg  and  pgWp, 
respectively.  F^  and  F2  represent  the  flux  terms  on  the  right-hand  sides  of 
the  appropriate  gas  and  solid  phase  momentum  equations,  respectively,  while 
and  represent  the  gas  and  solid  phase  source  terms.  S1^  and  S'^ 
represent  all  the  gas  and  particle  phase  source  terms,  except  for  the  drag. 
Finally,  0  =  (pg/4r2p)  fpg  from  the  drag  relation. 

Figure  28  shows  the  improvement  in  stability  caused  by  the  use  of  the 
predictor/corrector  strategy  on  the  gas  phase  z  -momentum  equation.  The 
momentum  is  plotted  for  a  single  node  versus  number  of  integrations.  Note 
that  the  kink  at  six  integrations  is  caused  by  the  particles  at  the  adjacent 
node  igniting  and  not  by  instability. 


Typical  3 -  Dimensional  Volume  Cell 


□  Axial  Velocity  (w)  Node 
V  Azimuthal  Velocity  (7)  Node 
O  Radial  Velocity  (TT)  Node 
*+  Pressure,  Temperature,  Density,  Porosity  Node 


Figure  16.  Single  Three-Dimensional  Volume  Cell 
Utilizing  the  Staggered  Grid 
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Figure  17.  Comparison  of  Conventional  and 
Staggered  Grids 
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Figure  19.  Determination  of  the  Velocity,  u  at  a 
Scalar  Node 
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Figure  24.  Determination  of  the  Value 
of  W  at  i  +  1/2 


Figure  25.  Prescribed  Radial  Velocity  Profile  for 
Holes  Larger  Than  One  Node 
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Figure  28.  Comparison  of  z-Momentum  Showing  Improvement  in 
Stability  Causeu  by  Predictor  Corrector  Strategy 
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SECTION  V 


RESULTS 


1.  INTRODUCTION 

The  numerical  scheme  described  in  the  previous  chapter  was  incorporated 
into  a  FORTRAN-V  computer  code.  In  this  chapter  some  initial  calculations  are 
presented  which  clearly  show  that  unsteady  multi-dimensional  flow  can  be 
modeled  by  this  code.  While  the  output  presented  here  is  still  of  a 
preliminary  nature,  it  nonetheless  shows  the  suitability  and  potential  of  our 
model  for  predicting  the  rapid  transients  encountered  in  this  type  of  work. 

As  was  shown  by  the  one-dimensional  model  in  Section  II,  a  totally 
confined  bed  of  small  energetic  particles  exhibits  rapid  localized  increases 
in  pressure,  temperature,  velocity,  and  density  (see  Figure  4).  Pressure  can 
increase  by  five  orders  of  magnitude  in  less  than  0.1  milliseconds.  One 
should  therefore  expect  difficulty  in  developing  a  three-dimensional  depiction 
capable  of  simulating  such  severe  transients.  Since  the  inclusion  of  multi¬ 
dimensional  effect#  and  the  partial  confinement  boundary  condition  greatly 
complicates  the  problem,  it  was  decided  to  defer  the  inclusion  of  the  mass 
loss  through  the  container  walls  (the  alternate  boundary  condition)  until  the 
ability  of  the  model  to  simulate  the  totally  confined  multi-dimensional  flame 
spreading  was  proven.  Paragraph  2  discusses  the  input  for  the  baseline  case 
in  which  a  totally  confined  bed  was  initiated  in  such  a  manner  that  the  flame 
spreading  would  be  two-dimensional,  while  Paragraph  3  presents  the  results  for 
this  case.  Paragraph  4  presents  the  currently  available  results  for  the  more 
complex  partially  confined  case. 
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2.  BASELINE  CASE  DISCUSSION 

All  the  test  cases  reported  in  this  section  modeled  granular  beds  of 
explosive  with  a  fixed  bed  radius  of  0.5  cm  and  lengths  ranging  from  3  cm  to  6 
cm.  Although  these  dimensions  are  obviously  smaller  than  those  of  actual 
explosive-containing  warheads,  they  are  sufficiently  large  for  us  to  study  the 
developing  transients  which  appear  during  the  early  portion  of  the  event. 

Also,  while  the  code  is  capable  of  modeling  three-dimensional  flame  spreading, 
only  two-dimensional  flame  spreading  results  will  be  shown  here  for  ease  of 
data  presentation.  (Contour  plotting  only  allows  three  parameters,  one 
dependent  and  two  independent  variables,  to  be  plotted.) 

Initiation  of  the  granulated  explosives  was  achieved  by  assuming,  as 
before,  that  a  prescribed  high-temperature  profile  exists  at  one  end  of  the 
bed  (see  Figure  29  ). 

The  temperature  at  several  prescribed  points  in  the  profile  was  assumed 
to  be  above  the  bulk  ignition  temperature  of  the  solid,  thus  allowing  the 
particles  to  ignite  at  time  t  =  0.  In  this  manner  the  user  can  control  the 
initial  spatial  domain  in  which  the  flame  spreading  will  take  place.  One- 
dimensional  flame  spreading  (no  variation  of  the  parameters  in  the  ?  - 
direction  or  I  -direction)  can  be  produced  by  assuming  a  profile  one  shown  in 
Figure  29,  where  the  gas  and  particle  temperatures  vary  only  in  the  z- 
direction.  Two-dimensional  flame  spreading  in  the  ?  and  z  -directions  can  be 
produced  by  assuming  gas  and  particle  initial  temperature  profiles  which  decay 
to  the  ambient  temperatures  in  both  the  r  and  the  z  -direction.  Finally,  full 
three-dimensional  flame  spreading  can  be  produced  by  assuming  an  initial 
profile  which  varies  in  all  (r-,  a-,  and  z-)  directions. 
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Input  parameters  for  the  two-dimensional  baseline  (no  vent)  case  are 


given  in  Table  4: 


TABLE  4:  INPUT  PARAMETERS 


Bed  Length,  L 
Bed  Radius,  R 
Bed  Angle,  o 

Particle  Radius,  r_ 

„  . 

Porosity,  $ 

Particle  Density,  p 

Ambient  Gas  Temperature,  T_ 

30 

Ambient  Particle  Temperature,  Tp 

Initial  Gas  Pressure,  Pfl  ° 

Particle  Ignition  Energy,  Eign 

Space  Increment,  u 

Space  Increment,  ar 

Angular  Increment,  ae 

Initial  Gas  Viscosity,  ua 

yo 

Prandtl  Number,  Pr 

Gas  Constant,  R' 

Specific  Heats,  C„  &  C„ 

9  P 

Particle  Chemical  Energy,  Echem 
Burning  Rate  Index,  n 
Burning  Rate  Coefficient,  a 
Gas  Equation  of  State  Constant, 
Number  of  Nodes 

Number  of  Variables  in  Storage 


3  cm  <  z  <  6  cm 
0,5  cm 


100  urn  <  rD  <  200  um 
0.30 
1.675  g/cm^ 
300  K 
300  K 
100  kPa 
4.6206  x  10®  ergs/g 
1  mm 


1.8  x  10  4  dyne-sec/cm^ 

0.7 

296.79  ergs/K  gmole 
1.51  x  10^  ergs/g  K 
5.74  x  1010  ergs/g 
1.0 
0.001 
4.0 

5  x  1  x  60 


3.  RESULTS  FOR  THE  TWO-DIMENSIONAL  SIMULATIONS 

Figures  30-46  show  the  pressure,  temperature,  radial  velocity  and  axial 
velocity  profiles  for  the  gas  at  various  times.  The  explosive  is  6  an  long, 
composed  of  200  ym  sized  particles.  The  grid  dimensions  were  5  by  1  by  60 
nodes,  and  the  bed  was  initiated  in  such  a  manner  that  the  flame  spreading  was  * 
two-dimensional  (in  the  ?-  and  z  -direction). 

Figures  30-33  show  the  rapid  increase  in  the  gas  pressure  at  the 
initiated  end  of  the  bed,  as  well  as  the  forward  motion  (along  the  z  -axis)  of 
the  pressure  front  as  time  increases.  (Note  that  the  scale  on  the  pressure 
axis  changes  on  each  figure.)  Particularly  interesting  is  the  peak  forming  at 
the  head  of  the  pressure  wave  at  t  =  46.3  usee.  This  is  indicative  of  a  shock 
forming  in  the  granular  explosive  bed. 

Calculations  were  not  possible  after  about  50  usee  (350  integrations)  due 
to  the  rapidly  developing  shock  described  above.  This  was  caused  by  the  fact 
that  we  are  using  a  constant  time  step,  rather  than  a  time  step  which  adjusts 
itself  to  the  rapidly  changing  bed  conditions.  As  the  shock  forms,  pressure 
gradients  across  the  front  of  the  wave  become  very  steep  and  much  smaller  time 
steps  are  needed  to  allow  the  differencing  scheme  to  adjust  to  these  strong 
gradients.  The  instability  was  caused,  then,  not  by  the  finite  differencing 
scheme,  but  by  the  simplification  of  assuming  a  constant  time  step. 

Most  interesting  of  all  is  to  observe  that  as  time  increases,  gradients 
in  the  radial  direction  vanish  entirely  and  the  flow,  which  was  initially  two- 
dimensional,  becomes  totally  one-dimensional  (with  gradients  occurring  only  in 
the  axial  direction).  This  is  exactly  what  the  experiments  of  8ernecker  and 
Price  Reference  9  have  shown  to  occur  when  totally  confined  explosives  are 
initiated  in  a  multi-dimensional  manner.  Figures  38-43  show  this  phenomena 
more  clearly. 
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Figures  34-37  show  the  gas  temperature  profiles  at  various  times  during 
the  burn.  One  can  clearly  see  the  gradients  in  the  radial  direction  during 
the  early  stages,  as  well  as  the  disappearance  of  these  gradients  during  the 
latter  stages. 

Figures  38-43  show  the  decay  of  the  radial  gradients  most  dramatically, 
since  the  disappearance  of  these  gradients  forces  the  radial  velocity  to  be 
zero.  In  Fig.  38  we  see  a  very  prominent  radial  velocity  profile  at  t  *  13.3 
usee.  The  radial  velocity  here  is  positive,  meaning  that  the  direction  of  the 
flow  is  outward  from  the  boundary  at  r  =  0.  Since  there  is  an  impermeable 
wall  at  r  =  R,  we  know  that  the  gas  must  eventually  reflect  off  of  this  wall 
and  reverse  its  direction. 

In  Figure  39  we  see  that  the  location  of  the  peak  radial  velocity  has 
moved  downstream  (in  the  z-direction)  and  the  velocity  behind  the  peak  is 
negative,  as  expected,  due  to  the  reflection  of  the  flow  in  that  region  from 
the  wall.  The  reverse  flow,  though  difficult  to  see  in  Figure  39  because  of 
the  manner  in  which  the  figure  is  plotted,  is  more  clearly  evidenced  by  the 
trough  behind  the  velocity  peak  shown  in  Figure  40.  In  Figure  41  we  again  see 
the  reverse  flow  trough,  noting  that  it  has  now  moved  farther  downstream. 

This  means  that  a  portion  of  the  flow  which  was  positive  in  Figure  40  (at  t  = 
34.47  ysecs)  has  reflected  off  the  wall  at  r  =  R  and  is  now  heading  back 
toward  the  center  of  the  bed  by  the  time  t  =  39.96  usee.  We  can  also  observe 
the  slight  peak  at  z  =  I  cm  where  a  trough  had  occurred  in  the  previous 
figure.  This  indicates  that  the  reverse  flow  from  Figure  40  has  again 

reflected,  this  time  from  the  boundary  at  r  =  0,  and  is  now  again  heading  in 

/ 

the  positive  ?  -direction.  One  should  also  note  the  overall  reduction  in  the 
magnitude  of  the  velocities,  indicating  that  the  flow  is  indeed  damping  out. 
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r.  w.  r~w. 


(Note  the  change  in  scale  on  the  radial  velocity  axis  as  compared  to  the 
previous  figures.)  Figures  42  and  43  show  this  behavior  continuing  until,  at 
time  t  =  46.36  usee,  there  is  practically  no  radial  velocity  in  the  bed. 

(Note  again  that  the  scale  on  the  radial  velocity  axis  changes  in  these 
figures.) 

Figures  44-47  show  the  developing  axial  velocity  profiles.  Again  we  note 
the  initially  two-dimensional  profile  decaying  to  a  one-dimensional  profile. 

In  Figure  48  the  location  of  the  flame  front  is  plotted  in  r,  z-space  at 
various  times  during  the  event.  We  can  note  that  initially  the  flame  front 
location  varies  in  both  the  r  and  z  -directions.  As  time  progresses,  the 
variation  of  the  flame  front  location  as  one  moves  in  the  ?  -direction  becomes 
less  and  less,  until  at  time  t  =  47.7  ysec  no  variation  in  the  radial 
direction  is  predicted.  Clearly,  the  flame  front  has  become  one-dimensional. 

These  results  illustrate  two  key  facts.  First,  the  fact  that  the 
numerical  scheme  described  in  Section  IV  was  able  to  handle  the  rapid  flow 
transients  indicates  its  suitability  for  use  in  analyzing  dynamic  two-phase, 
multi-dimensional  flows.  Secondly,  the  results  also  show  that  the  model 
developed  in  Section  III  predicts  behavior  which  is  at  least  in  qualitative 
agreement  with  data  from  experiments  involving  accelerating  reaction  fronts  in 
granular  beds  of  explosive  solids. 

One  may  ask  at  this  point  if  a  multi-dimensional  flame  spreading  model  is 
only  of  limited  use,  since,  in  the  presence  of  total  confinement,  the  flame 
spreading  will  eventually  asymptote  to  a  one-dimensional  process.  One  must 
realize  that  while  the  predictions  presented  in  Figures  30-46  show  that  multi¬ 
dimensional  effects  last  only  50  usee  for  our  example  baseline  case,  such  an 
amount  of  time  represents  a  significant  percentage  of  the  total  duration  of 
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typical  DDT  events,  which  may  take  only  100-200  ysec.  Hence,  this  transient 
two-dimensional  stage  will  indeed  have  a  significant  effect  on  a  prediction  of 
the  detonation  run-up  length. 

Figures  49-51  show  the  gas  pressure,  axial  velocity,  and  flame  front 
profiles  in  a  bed  which  was  initiated  in  a  purely  one-dimensional  manner,  and 
compare  them  to  the  results  from  the  same  bed  initiated  in  a  two-dimensional 
manner.  In  both  of  these  case  the  initial  particle  size  was  fixed  at  a  100  ym 
radius.  Figure  49  shows  the  gas  pressure  profile  at  t  =  19.5  usee  for  the 
case  in  which  the  initiation  was  one-dimensional  compared  with  the  two- 
dimensional  case  at  r  =  R/2,  again  at  t  =  19.5  ysec.  One  can  see  the  obvious 
difference  in  the  two  profiles.  Figures  50  and  51  provide  similar  comparisons 
of  the  axial  velocity  and  the  flame  front  locus,  respectively.  Again,  the 
profiles  are  markedly  different  for  the  two  different  types  of  initiation. 
Therefore,  we  can  conclude  that  while  the  multi-dimensional  effects  eventually 
decay  they  still  have  a  significant  effect  on  the  results. 

4.  PARTIAL  CONFINEMENT  CASE 

Once  the  capabilities  of  the  code  had  been  shown  for  the  baseline  (no 
vent)  case,  calculation  was  made  to  test  the  model's  ability  to  simulate  flow 
with  partial  confinement.  Again  a  3  cm  long  bed  of  particles  with  an  initial 
radius  of  100  ym  was  modeled.  A  crack  3  mm  (3  nodes)  in  width  was  assumed  to 
exist  around  the  entire  circumference  of  the  bomb.  This  latter  condition  was 
imposed  to  cause  the  flow  inside  the  domain  to  effectively  remain  two- 
dimensional.  The  crack  was  located  a  distance  of  t  =1.55  cm  from  the 
initiated  end  (z=0).  A  uniform  profile  was  used  to  assure  that  the  initial  * 
flame  spreading  one-dimensional,  in  order  to  make  the  multi-dimensional 
effects  caused  by  the  hole  located  at  a  fixed  radius  to  be  more  apparent. 
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Figure  52  shows  ,,ie  predicted  gas  pressure  profile  at  t=21.48  usees.  One 
can  see  that  behind  the  hole  the  profile  is  indeed  one-dimensional.  The  small 
“indentation"  in  the  profile  at  z=1.55  cm  indicates  the  location  of  the 
hole.  The  flow  at  this  point  is  still  largely  one-dimensional  because  the 
pressure  front  has  only  just  reached  the  hole. 

Figure  53  shows  the  gas  pressure  profile  at  a  time  1.5  visecs  later.  By 
then  the  effect  of  the  hole  is  very  apparent.  The  pressure  is  severely 
diminshed  at  the  crack  itself,  while  downstream  of  the  hole  (i.e.,  in  the 
positive  7  -direction)  one  can  see  strong  gradients  in  the  r  -direction, 
indicating  that  the  flow  is  multi-dimensional.  It  is  especially  interesting 
to  note  that  downstream  of  the  hole  the  pressure  is  higher  at  the  wall  than  at 
the  centerline.  This  is  caused  by  the  fact  that  downstream  of  the  hole,  the 
radial  velocity  at  r=R  (i.e.  the  solid  wall)  must  again  be  equal  to  zero.  The 
radial  velocity  components,  induced  by  the  gas  flowing  out  of  the  hole,  must 
decelerate  again  to  zero  at  the  wall,  thus  increasing  the  pressure. 

Figure  54  shows  a  plot  of  the  radial  velocity  profile  shortly  after  the 
pressure  wave  has  reached  the  hole.  One  can  see  that,  due  to  the  one- 
dimensional  initiation,  the  radial  velocity  is  zero  everywhere  except  at  the 
hole.  In  Figure  55,  the  pressure  signals  from  the  hole  have  propagated  down 
to  r=0  and  upstream  of  the  hole,  inducing  radial  velocity  components  in  these 
regions. 

In  Figure  56  one  can  see  that  the  velocity  profile  is  even  more 
prominent,  with  a  dimple  in  the  profile  upstream  of  the  crack.  This  Indicates 
that  the  radial  flow  induced  by  the  hole  has  reflected  from  the  solid  wall  at 
r=R  and  is  flowing  back  toward  the  centerline. 
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For  this  partially  confined  case,  the  calculations  were  continued  until 
tne  flame  front  reached  the  end  of  the  bed  at  t=23.2  usee,  at  which  point  the 
calculations  were  automatically  shut  off.  This  is  somewhat  significant  in 
that  when  the  code  was  run  for  the  identical  input  conditions,  except  with 
total  confinement,  the  formation  of  a  shock  wave  prevented  further 
compulsions  after  a  time  of  19.0  nsec.  In  other  words,  the  existence  of  the 
hole  sufficiently  weakened  the  formation  of  a  I  -direction  shock  so  that  the 
scheme  remained  stable  for  the  same  spatial  griding. 

The  results  presented  in  this  section  clearly  indicate  that  the  model  and 
code  developed  in  this  study  has  the  capabilities  needed  for  calculating 
partially  confined,  unsteady  multi-dimensional  flows.  The  potential  for 
carrying  out  calculations  for  (more  realistically)  larger  domains,  with 
different  vent  openings  and  for  explosives  with  different  granulation  is  very 
clear. 

5.  SUMMARY  ANO  RECOMMENDATIONS 

Two  models  were  developed  to  predict  the  behavior  of  unsteady  two-phase 
reactive  flows  in  a  partially  confined  region.  A  pseudo  two-dimensional  model 
was  applied  for  a  wide  variety  of  cases  and  the  results  indicated  that  partial 
confinement  does  have  a  significant  quenching  effect  on  a  forming  detonation 
wave. 

Since  the  pseudo  two-dimensional  formulation  required  many  major 
simplifying  assumptions,  it  was  decided  that  a  three-dimensional  unsteady 
model  should  be  developed  to  provide  a  more  realistic  fluid  mechanics 
simulation.  The  conservation  equations  were  derived  in  cylindrical 
coordinates  and  a  new  numerical  scheme  was  devised  to  integrate  these 
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equations.  At  the  time  of  this  publication,  preliminary  calculations  of  this 
model  indicates  that  it  will  provide  stable,  physically  meaningful  results 
which  are  consistent  with  the  very  limited  experimental  evidence  available. 

Time  constraints  on  the  publication  of  this  analysis  prevents  a  larger 
number  of  calculations  from  being  performed  with  this  multi-dimensional  model. 

To  conclude,  the  following  proposed  improvements  are  recommended  for 
future  work. 

(1)  Derivation  of  an  Optimal  Time  Step: 

In  order  to  be  able  to  predict  the  formation  of  a  steady-state 
detonation  wave  within  the  fragmented  explosive  bed,  stability  of  the 
numerical  differencing  model  must  be  maintained  in  spite  of  the 
presence  of  a  shock  wave  within  the  computational  domain.  To  achieve 
this,  the  time  step  must  be  sensitive  to  the  instantaneous  conditions 
inside  the  bed.  A  Von  Neumann-type  stability  analysis  (Reference  16) 
may  be  necessary  to  determine  the  formula  for  the  optimum  time  step. 

t 

(2)  Incorporation  of  Improved  Constitutive  Relations: 

As  noted  in  Appendix  A,  the  applicability  of  the  data  required 
for  the  constitutive  relations  under  the  wide  range  of  conditions 
imposed  by  this  problem  is  questionable.  The  incorporation  of  the 
new  gas  equation  of  state  be^ng  developed  by  Wang  and  Krier  Reference 
22  and  the  Mie-Grunniesen  solid  phase  porosity-pressure  relation 
should  help  improve  the  accuracy  of  the  model. 

(3)  Devising  Alternate  Forms  of  Initiation: 

The  use  of  a  prescribed  temperature  to  ignite  the  particles, 
while  proving  satisfactory  for  initiating  the  bed,  is  not  always 


realistic.  Explosive  warheads  which  are  to  be  modeled  are  often 
initiated  by  means  of  a  shock  from  an  exploding  fuse.  Clearly,  at 
time  t=0,  not  only  temperature  but  high  local  pressure  should  be 
modeled. 

(4)  Application  of  New  Generation  Computer  Technology: 

The  CDC  CYBER- 175  computer  currently  in  use  at  the  University  of 
Illinois  at  Urbana-Champaign  was  employed  to  run  the  code  shown  in 
Appendix  D.  While  this  is  one  of  the  fastest  and  largest  mainframe 
computers  on  the  market  it  is  not  totally  adequate  for  all  the 
aspects  of  this  type  of  work.  Though  we  have  been  able  to 
demonstrate  the  capabilities  of  the  code  for  predicting  multi¬ 
dimensional  flame  spreading  in  small  domains  using  the  CYBER-175,  the 
storage  and  processing  requirements  for  applying  this  code  to  larger 
domains  (with  full  three-dimensional  flame  spreading)  would 
necessitate  the  use  of  an  even  more  powerful  computer. 

The  code  requires  that  the  values  of  ten  primary  variables  be 
stored  at  two  time  levels,  as  well  as  the  values  of  twelve  secondary 
variables  at  one  time  level,  at  each  node.  The  core  memory  of  the 
CYBER-175  allows  the  user  131,000  words  of  storage.  Leaving  enough 
memory  for  the  code  itself,  this  restricts  the  maximum  computational 
domain  to  3500  nodes.  Given  the  size  of  the  space  increments  needed 
to  assure  proper  resolution  (approximately  1  x  10"3m),  this  limits  us 
to  calculating  cases  representing  relatively  small  domains.  Of 
course,  infinite  storage  is  available  by  using  off-line  memory,  but 
only  at  an  increased  cost  to  the  user. 
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This  problem  may  be  solved  by  the  utilization  of  the  new 
generation  super-computers,  such  as  the  CRAY-1.  The  increase  in  core 
memory  (1-4  million  words  for  the  CRAY-1)  offered  by  such  a  computer 


would  allow  us  to  perform  calculations  for  much  larger  domains. 
Furthermore,  the  increased  computational  speed  of  these  computers 
will  help  keep  the  computational  costs  to  a  minimum. 
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Figure  30.  Gas  Pressure  Profil 
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Figure  32.  Gas  Pressure  Profile  at  t  =  39.96  usee 
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Figure  36.  Gas  Temperature  Profile  at  t  =  39.96  ysec 
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Figure  39.  Gas  Radial  Velocity  Profile  at  t  -  28.06  usee 
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Gas  Radial  Velocity  Profile  at  t  =  39.96  nsec 
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Figure  43.  Gas  Radial  Velocity  Profile  at  t  =  46.36  ysec 
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Figure  45.  Gas  Axial  Velocity  Profile  at  t  =  34.47  usee 
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Figure  46.  Gas  Axial  Velocity  Profile  at  t  =  39.96  psec 


Figure  48.  Flame  Front  Locus  Showing  Transition  from  Two-Dimensional  Flame 
Spreading  to  One-Dimensional  Flame  Spreading 
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Figure  51.  Comparison  of  the  Flame  Front  Locus  for  One- 
and  Two-Dimensional  Initiations 
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Figure  52.  Gas  Pressure  Profile  at  t  -  21.48  usee 
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Figure  53.  Gas  Pressure  Profile  at  t  =  22.97  ysec 
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APPENDIX  A 


CONSTITUTIVE  RELATIONS 


I 


i 


To  totally  describe  the  process  occurring  in  the  combusting  explosive 
media,  we  must  be  able  to  find  the  instantaneous  values  of  pg,  pp,  ug,  up,  Tg, 
Tp,  Pg,  Pp  and  $.  The  six  conservation  equations  allow  us  to  solve  for  six  of 
these  variables.  In  addition  to  these  equations,  three  constitutive  relations 
are  needed  for  closure: 

(i)  An  Equation  of  State  for  the  Gaseous  Products  Phase: 

In  this  work  we  utilized  a  non-ideal  equation  of  state  for  hard  spheres 
suggested  by  S.  J.  Jacobs  (Reference  (231).  This  equation  of  state  was 
written  as: 

Pg  ■  »gR'T,  0  *  Vg>  <A_1) 

Ideally,  the  coeficient  b^  should  be  a  function  of  the  density.  However,  in 
this  analysis,  b^  was  considered  a  constant. 

While  offering  good  correlation  with  experimental  values  at  relatively 
low  pressures,  there  is  some  question  about  the  accuracy  of  this  equation  of 
state  at  the  extremely  high  pressures  (10-20  GPa)  which  can  be  achieved  in  a 
bed  of  combusting  particles.  Work  is  currently  being  performed  by  Krier  and 
Wang  (Reference  [22])  to  develop  an  equation  of  state  which  is  va14d  at  these 
extreme  pressures. 

(ii)  An  Equation  of  State  for  the  Solid  Particle  Phase: 


Here,  the  Tait  equation  of  state  was  used  to  predict  the  density  of  the 
solid  particles: 


(iii)  A  Relation  for  Predicting  the  Porosity,  $: 

Since  no  relationship  for  the  porosity  of  a  packed  bed  of  particles  which 
is  reliable  under  the  range  of  conditions  imposed  by  this  problem  is  known, 
this  condition  was  replaced  by  an  equilibrium  condition  which  states  that  at 
each  instant  of  time,  the  particle  pressure  is  equal  to  the  gas  pressure.  The 
porosity  is  then  predicted  for  each  new  time  step  by  taking  the  value  of  o ^ 
predicted  by  Equation  (2)  and  dividing  it  by  the  solid  phase  density  from  the 
previous  time  step.  While  this  procedure  does  introduce  some  error  into  the 
prediction  of  *,  this  error  will  be  quite  small  since  the  time  steps  used  in 
this  model  are  quite  small  and  the  density  of  the  solid  phase  will  not  change 
much  between  adjacent  time  steps. 

in  addition  to  these  equations,  relations  are  needed  to  solve  for  the 
heat  transfer,  drag,  and  burning  rate. 

( i v )  A  Relation  for  the  Interphase  Heat  Transfer  Term,  Q  : 

The  relations  used  to  find  the  heat  transfer  between  the  hot  gases  and  the 
particles  were: 


(A— 3) 


Q  =  h  (T  -  T  }  .3ikiL 
P9^  9  P1  rp 

h  *  -4  (1  +  0.2  Re  0,7  Pr0'33) 
P9  *■  r  ' 


(A-4) 


where  the  gas  thermal  conductivity  is  defined  as: 


kg  *  “g<R'+  C.g>/Pr 


(A-5) 


and  the  interphase  viscosity  is  given  by: 


0  6  5 

u  =  U  (T  /T  )  * 

g  gn  g  gft 


( A-6) 


The  Reynolds  Number  used  in  this  work  is  given  by  Wallis  [241  for  two-phase 
flow  as: 


Rer  ■  2  V‘I(up  "  Ug)|/Bg 


(A-7a) 


For  the  three  dimensional  model  described  in  Section  II,  a  variation  of  this 
form  of  the  Reynolds  Number  was  used: 

Rer  =  2rp0l|(up  +  vp  +  wp  )  -  (ug  +  vp  +  wp  )  |/ug 

(A-7b) 

(v)  A  Relation  for  the  Gas-Particle  Orag,  D: 

The  relation  used  for  the  interphase  drag  in  this  work  are  given  as: 


0  *  »gl(Ug  -  >1'  fpg 


fpg  ■  t1? !  ‘  (226  *  51^1  °-*7) 


(A— 8) 

(A-9) 


This  relation  was  derived  experimentally  by  Kuo  and  Nydegger  (Reference  [25]) 
for  steady  flow  in  a  packed  bed  of  spheres  at  constant  porosity  for  Reynold's 
number  ranging  from  460  to  14,600.  The  problem  being  considered  here  is 
unsteady  and  has  Reynold's  numbers  which  greatly  exceed  this  range.  However, 
at  the  present  time  there  does  not  appear  to  be  any  relationship  in  use  which 
is  any  more  rel iable. 

(vi)  A  Relation  for  the  Gas  Production  Rate,  r: 

r  =  ^  (1  -  •)  or  (A- 10) 

P  v 

where  rp  is  the  instantaneous  particle  radius  and  r  is  the  surface  burning 
rate  given  by: 

r  »  a  Pgn  (A— 11) 

While  some  of  the  relations  used  in  this  work  were  applied  under 
conditions  which  lay  outside  their  stated  domain,  this  should  not  detract 
anything  from  the  reliability  of  the  results  obtained  herein.  Rather,  the 
development  of  constitutive  relations  which  are  tailored  to  the  extreme  range 
of  conditions  encountered  in  this  type  of  work  will  merely  help  to  fine  tune 
the  previous  results. 


APPENDIX  B 


THE  SPEED  OF  SOUND  IN  A  NON-IDEAL  GAS 

In  Section  II  we  formulated  a  partial  confinement  model  which  assumed 
that  the  gas  escaped  from  each  damaged  control  volume  at  the  choked 
velocity.  We  must  then  determine  the  speed  of  sound  in  the  non-ideal  gas 
which  is  produced  in  the  .burning  particles.  Recall  that  our  non-ideal 
equation  of  state  had  the  form: 

P  =  pR'T(l  +  blP)  (B-l) 

Recall  also  that  the  First  Law  gives  us: 

de  =  «Q  -  pdv  (B-Z) 


If  we  use 


de  =  CydT 


and  the  Second  Law 


ds  =  f 


we  can  re-write  (B.2)  as 

i 

« 

CydT  =  Tds  -  pdv 


(B-3)* 


(B— 4) 


(8-5) 
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*  -  (If),  dT  +  (ff)T 


and  applying  the  First  and  Second  Laws  of  Thermodynamics  and  using  Maxwell's 
relations. 


i 


1 
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APPENDIX  C 

DERIVATION  OF  THE  CONSERVATION  EQUATIONS 

In  Section  III,  the  ten  conservation  equations  for  the  unsteady,  three 
dimensional  separated  flow  model  were  presented  in  their  finished  form.  In 
this  appendix,  the  key  steps  in  deriving  these  equations  will  be  presented. 


Gas  Phase  Conti nui tv 


Consider  a  three  dimensional  cylindrical  control  volume  element  with 


sides  of  finite  length  dr,  da,  and  dz.  Mass  is  flowing  through  the  volume  in 
the  r,  8,  and  z  directions.  Performing  a  mass  balance  on  this  control  volume: 


Mass  increase  in  C.V. 
in  time  t 


net  mass 


mass  generated 


flux  through  +  in  C.V.  in  time  t 
C.V.  in  time  t 


'  <A>r  -  <">r.dr  +  <*>8  -  K+de*  rcvT 

(C-l) 

Here,  Vg  =  volume  of  gas  phase  and  Vy  =  total  volume  of  control  cell.  Using  a 
Taylor  Series  to  expand  the  flux  terms  about  r+dr,  o+do,  and  z+dz,  we  get 


3(p„VJ  2 

—ft3"  =  V  t"V+  “  dr  *  “V  dr  + - > 


.  2  . 

3  m  3  m  2 

+  V  (V  — rr  da  +  — ^  da  +  -  - 
80  39  30Z 

.  2  . 
am  a  m  2 

+  -  {m  +  dz  + - 1  dz  +  - 

3Z 
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or 


3t 


2. 

am,.  a  m  2 

b? dr  ♦  -4  <*• 

ar 


+  -  -  - 


am. 


■39 


2  . 
3  m 


de  + 


ae 


4  de2  + 


2  . 

am  am, 

-  (it  dz  +  — 4 

32 


-  } 

-  -  } 
-  -  } 


+  rV 


T 


(C-2) 


For  a  cylindrical  control  volume  Vg=r,drdodz$ .  Since  the  volume  of  the  control 
cell  is  invariant  with  time,  it  may  be  pulled  out  of  the  time  derivative  on 
the  left  hand  side.  Oividing  by  Vg,  we  get 


a(cq$) 

3t 


2. 

,  3«L  3  m 

L  bS *  — ?  <*-  + - i 


rdedz  1  3r 

1 


ar 


rdrdz  l39 
1 


.  2 . 
am„  3  rn 

(_1  +  9  d9  + - } 


am. 


rdrde  laz 


{- 


30 
2  . 

3  m. 


3Z 


dz  +  }  +  r 


(1) 

(2) 

(C-3) 


How  consider  the  first  group  of  terms  (1)  on  the  right-hand  side  (R.H.S.)  of 
Equation  (C-3): 


1 

rdedz 


dr  +  - - } 


Now: 
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m 

A 


r 

r 


p  A  u 
g  r  g 

rdedz$ 


Therefore,  we  get: 

2 

-  7353?  <1?  (0  rdsd2«u  )  «■  -iy  (»  rdedz*u  )  dr  * - } 

3  r 

Since  de  and  d2  are  not  functions  of  r,  they  may  be  pulled  out  of  the 
derivative.  Therefore,  we  get: 

-  r  <77  (o'rug)  *  <“irug)  dr  + - -  ) 


By  defining 


and 


mn  =  pA  v„ 
e  g  6  g 


An  =  drdz$ 


0 


m_ 


pgAzWg 


A  ^  *  rdrd8<* 


and  performing  the  same  analysis  on  groups  (2)  and  (3),  we  finally  get 


a(pi) 


-  7  <77  <r»>ug)  +  A  <r°iug)  dr  + 

y  a  r  y 


r  lar 


-  -  -}  (1) 


at 


2 


+  r 


PlVg)  d6  + - }  (2) 

Vg)  dz  +  -  -  -  -}  (3) 

(C— 4) 


If  we  now  take  lim  dr  -  0,  all  of  the  higher  order  terms  in  group  (1)  approach 
zero.  Similarly,  if  we  take  lim  de  -  0  and  lim  dz  -  o,  all  the  higher  order 
terms  in  groups  (2)  and  (3)  go  to  zero.  Finally  then,  we  are  left  with: 


ap  t 

~Tt 


far  <ro  'ugj  ’  7  25  (°.vg>  'SI  <’.wg>  + 


(C-5) 


Gas  Phase  Momentum  Equations 

Since  velocity  is  a  vector  function,  the  momentum  equations  must  be 
derived  in  vector  form.  Consider  again  a  cylindrical  control  volume. 
Performing  a  momentum  balance  on  this  control  cell: 

Momentum  increase  net  momentum  net  pressure  interphase 
in  C.V.  in  time  t  =  flux  through  +  stress  on  C.V.  +  drag  and 

C.V.  in  time  t  faces  momentum  generated 

by  combusting  particles 

*  Vg5)  ■  (IW  <Ar»W  (tfU. 

<PgAr>r-  <Pg\W 
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where 


+  (P  A  )  -  (P  A  )  ,  +  (P  A  )  -  (P  A  ) 

1  g  a'b  v  g  e'9+de  1  9  z'z  v  g  z;z+dz 


D  V  +  r  7  VT 
g  F 


(C-6 


7=ur+v5+w2 

g  g  g 


D  =  D  +  D  +  0 
r  9  z 


Again  performing  a  Taylor  Series  Expansion  as  before,  we  get:. 


-i  SV1  • 


fc?  (■-»)  dr  +  (i  f)  dr  +  -  -  }  (4) 
a  r 


^rf?  rde  +  ” — 2  (%^)(rd0)^+  *  -  1 


(rde) 


-  {gf  (M)dz  +  (mzV)  dz  +  -  -  }  (6) 

9Z 


-  {~  (PA)rdr  +  -ij  (PA)  dr  +  -  -  }  (7) 

ar 


Ifi?  <pSVde  + 


5  (PS)  Jrde)  ♦  -  -  } 


(rse) 


(PA)  dz  +-^7  (PA)  dz  +  -  -  }  (9) 

3Z 


-  D  V  +  r  V  9 

g  g 


(C-7 


Let  us  consider  each  term  separately. 


Term  4: 


-+ 

a(nrV) 

— —  =  7-  [m^(u„r  +  v  e 
3 r  ar  1  r'  g  g 


=  — ^  [m  ur  +mv$+mwz 
ar  1  r  g  r  q  r  q 


(C-8) 


Nate  that  —  =  =  0  .  Therefore, 

3 1  01  dr 


4-  (m  V) 
ar  x  r  ' 


-  7-  (m  u  )  ?  +  -4  (m  v  )  e  +  -d-  (m  w  )  2 
ar  '  r  a'  ar  v  r  n'  ar  v  r  qJ 


Now, 


mr  =  PgUgrdadz$ 

Therefore, 

aF  (li,rV)  =  d9dz  Ir  <°irug2)  ?  +  d0dz  If  (pxugvgr)  9 

+  dedz  (piugwgr>  z  (c~9) 

or,  ~  (ihrV)  =  dedz  (Plug2r)  ?  +  ^  (piugvgr>  5 

+  ar  (°!ugwgr)  z  1  (4*>  (C'l°) 

Note:  the  higher  order  terms  will  not  be  expanded,  since  they  will  disappear 
as  before.  Going  on  to  term  5,  we  proceed  in  a  like  manner,  recalling, 
however,  that  ||=S,||=-f,-|4*0.  Therefore, 


Therefore, 


A  =  rdedzs  r 
r 


(pgAr)dr  =  (rdedzPg$)  r  =  drdedz  (P,r)  f 

There  is  another  component  of  the  pressure  stress  in  the  r  -direction 
which  is  due  to  the  curvilinear  nature  of  our  coordinate  system.  Consider  the 
r,9  projection  of  our  control  volume  shown  in  Figure  57.  As  can  be  seen  from 
the  drawing,  there  is  a  component  of  the  pressure  force  in  the  $  -face  which 
acts  in  the  r  -direction.  The  magnitude  of  this  force  is: 

P 

PgAe  Sin  (^)  =  PgAe  <^f)  =  “§  drd8dz*  (C-14) 

Since  an  equal  force  is  acting  on  the  opposite  §  -face,  the  total  force  is: 

fy  =  -Pgdrdedz$r  =  -Pjdrdedzr  (C- 15) 

Therefore,  term  7  becomes: 

~  (PgAr)r  =  drdedz  (Ptr)  -  Pt)  f 

3P! 

=  drdedz  (Pt  +  r  — —  P,)  r 


(Pg^r)r  =  rdrdedz  (P,)|  r  (7*)  (C-16) 


Proceeding  to  term  8 

rfe  <PgBa>erde  '  lPgdrdz»l  rd8  5 

»  rdrdsdz  (-j|  (P,)|  0  (C-17) 

Again  there  i$  component  of  the  stress  on  the  f  -face  in  the  a  direction: 

-  PgSin  $•)  Ar  +  dr  ♦  PgSin  (^|)  Ar  -  F  ...  (C-18) 

t  r-5  "  Pj^f  d0dz*  <r  '  (r+dr)) 

F  -  t  *  -  ijP  drdeZdzo 
r-e  g 

Therefore,  term  8  becomes: 

(PgAp0rde  a  drdzde  Ijf  (p.)  -  \  p.del  5  <8*>  (c-19> 

Similarly,  term  9  becomes: 

af  (Pg*2)z  dz  =  rdrdadz  £  (P,)  (9*)  (C-20) 

Now,  we  plug  these  starred  terms  back  into  Equation  C.4  and  divide  out  the 
volume  of  the  control  cell  and  take  the  limits  as  dr,  da,  and  dz  -»  0.  After 
re-arranging  the  terms,  we  finally  get: 
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Gas  Phase  r  Direction  Momentum 


it  <°iuq>  ‘  '  r  ir  '  7  if  <».yg>  '  it  (o>Vg) 


r  ar  1  g 

2 


o,v 


+  — 3 - »  (P.)  -  D  + 

r  ar  v  i'  r 


ru_ 


(C-21) 


Gas  Phase  8  -Direction  Momentum 


it  (»'V  =  -  r  if  ("-“aV'  -  7 if  <»>y  -  it  (“.Vs1 


-  .  1  _»  (P  )  .  o  +rv 

r  r  39  '  9  t 


(C-22) 


Gas  Phase  z  -Direction  Momentum 


it  (“iwg)  ‘ 


7  if  (».Vgr)  ’fit  (,,'vgV 


it  (“."o’)  'it  (pi>  -  0,  ♦  r“r 


(C-23) 


Gas  Phase  Energy  Equation 

Once  again  we  perform  an  energy  balance  on  our  control  volume: 

increase  in  energy 
in  control  volume 
over  time  t 


net  energy  flux  through 


control  volume  +  work  done 
in  moving  gas  through  C.V. 


heat  transfer,  drag  work, 
chemical  energy,  and  kinetic 
energy  of  the  combusted  particles 


3  ( o  E  V  )  P  P 

,  [A  (E  *_£)]  .  (A  (E  .+  -a)i  . 
at  1  r'  gt  og  Jr  1  rv  gt  og  Mr+dr 

♦  ivvjiu 


-  QVt  -  Du  VT  -  D  v  V,-  D  w  VT 
I  r  p  i  apt  z  p  i 

2  7  2 

* r  <4™  *  Y 4  Y + ^f>  vr  <c-m> 


Proceeding  in  an  identical  manner  as  before,  we  finally  get: 


1  3  P- 


Tt  <».V  ■  -  7T?  I».V(V^)I  -  r  ai-  lp',g(V^)l 


~  Id,w  (E  «+  -9-)]  -  Q  -  0  u  -  D  v 
az  u 1  g  '  gt  p„  w  r  p  e  p 


2  2  2 

"  °zwp  +  r  (EJhe«  +  -|"  +  4-  +  Hr)  (C'Z5) 
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p  =*  c*fe.c\  r\or  ty\o.\ 
uw\\  Vtc^-or 


Figure  57.  Control  Volume  Showing  Component  0-Stres 
Acting  in  the  r-Di recti  on 
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